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solving  by  the  Newton-Raphson  technique.  The  results  yielded  nonlinear  limit 
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increased,  the  flutter  speed  decreased,  while  the  flutter  frequency  increased 
toward  the  torsional  natural  frequency.  The  analytic  results  compared  favorably 
with  previous  experimental  works  by  Dunn.  The  current  nonlinear  analysis 
procedure  seems  an  effective  technique  for  analyzing  stall  flutter  phenomena. 
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Warren  C.  Chen  entitled,  "A  Formulation  of  Nonlinear  Limit  Cycle 
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Dugimdji,  the  Principal  Investigator,  and  the  supporting  laboratory  staff. 
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CHAPTER  1 


INTRODUCTION 


Stall  flutter  deals  with  the  self-excited  oscillations  of  a  wing  in  a  separated 
or  stalled  flow.  The  classical  flutter  theory  is  traditionally  based  on  small 
amplitude,  smooth,  linear  potential  flow  and  is  well  understood  theoretically  [see 
Refs.  1  to  3]  The  nonlinear  behavior  of  the  large  amplitude  stall  flutter  motions, 
however,  involves  flow  separation  and  highly  nonlinear  aerodynamic  phenomena. 
This  makes  it  difficult  to  analyze  analytically,  and  one  must  generally  resort  to 
either  computational  methods  or  experimental  semi-empirical  methods  to 
characterize  such  behavior  and  predict  its  occurrence. 

Problems  of  stall  flutter  arise  in  connection  with  wings  at  high  angle  of 
attack.  If  the  wing  is  near  the  stall  region,  a  nonlinear  stall  flutter  limit  cycle 
oscillation  may  occur.  This  nonlinear  flutter  phenomenon  may  take  place  at  a 
lower  velocity  than  linear  theory  would  suggest.  Since  some  current  aircraft  are 
capable  of  maneuvering  at  high  angle  of  attack,  it  is  of  interest  to  explore  this 
nonlinear  stall  flutter  behavior. 

Modeling  of  the  dynamic  stall  phenomenon  has  been  a  primary  concern  in 
the  study  of  aeroelasticity.  Numerous  research  in  this  subject  were  done  for  over 
two  decades.  There  were  two  main  approaches,  theoretical  and  semi-empirical, 
which  bases  on  e}q)erimental  data. 


13 


14 


Discrete  potential  vortex  method  [Refs.  4  to  6],  a  theoretical  approach, 
ignores  the  viscous  terms  in  the  fundamental  equations  and  assumes  potential 
flow  without  the  boundary  layer.  The  zonal  methods  [Refs.  7  to  9],  also  a 
common  theoretical  approach,  under  certain  assumptions,  models  the  viscous, 
non-viscous,  and  transition  regions  of  the  flow  separately.  These  theoretical 
models  are  extremely  computational  intensive  and  are  limited  by  the 
approximations  of  their  formulation. 

The  semi-empirical  methods  attempt  to  use  static  data  with  corrections  to 
model  the  dynamic  stall  event.  The  method  only  models  the  gross  aspects  of  the 
phenomenon.  This  is  advantageous  because  the  static  data  already  takes  into 
account  some  of  the  aerodynamic  parameters  such  as,  the  effects  of  Reynold's 
number  and  airfoil  shapes.  There  were  numerous  researches  on  the  semi- 
empirical  analysis  [Refs.  10  to  21].  The  semi-empirical  method  is  generally  not 
computationally  intensive,  and  is  suitable  for  routine  aeroelastic  analysis. 
Extensive  work  on  creating  databases  of  static  data  was  done  by  McAlister,  Carr, 
&  McCroskey  [Ref  22]  for  the  NACA  0012  airfoil.  These  work  was  further 
extended  by  McAlister,  Pucci,  McCroskey,  &  Carr  [Refs.  23  and  24]  to  include  a 
wider  range  in  the  variable  parameters. 

The  specific  objectives  of  the  current  investigation  are  to  explore 
analytically  the  roles  of  nonlinear  aerodynamics  in  high  angle-of-attack  stall 
flutter  of  aircraft  wings,  while  attempting  to  develop  a  simple  nonlinear  method  of 
analysis  that  is  not  computationally  intensive  by  modifying  linear  theory. 

This  research  project  is  a  part  of  series  of  studies  at  the  Technology 
Laboratory  for  Advance  Composites  (TELAC)  at  M.I.T.  in  a  continuing  effort  to 
investigate  the  aeroelastic  flutter  behavior  of  aeroelastically  tailored  composite 
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aircraft  wings.  The  results  of  current  investigation  were  compared  with  previous 
experimental  work  done  at  TELAC  by  Dunn  [Ref.  14]. 

In  this  report,  Chapter  2  contains  the  description  of  basic  flutter  phenomenon 
and  reviews  the  linear  solution  methods  common  to  aeroelastic  analysis.  In 
addition,  description  of  the  preliminary  work  in  the  current  investigation  using 
linear  theory  solution  in  preparation  to  approach  the  stall  flutter  problem  are  also 
included. 

Chapter  3  describes  the  nonlinear  theory  which  the  current  investigation  is 
based  on.  Analytically,  this  chapter  seeks  to  expand  on  and  improve  the  efforts  of 
the  previous  investigations.  This  chapter  describes  the  simplification  of  modeling 
nonlinearity,  particularly  the  aerodynamic  formulations,  over  the  previous 
investigations. 

Chapter  4  details  the  results  of  the  theoretical  investigations,  comparing 
results  against  previous  works,  while  Chapter  5  gives  concluding  remarks  on  the 
significant  contributions  of  the  current  investigation  and  recommendations  for 
future  work. 


CHAPTER! 


REVIEW  OF  LINEAR  THEORY 


2.1  Flutter 

Aeroelasticity  deals  with  the  interaction  between  aerodynamic,  elastic  and 
inertial  forces.  The  aerodynamic  forces  are  external  forces  which  arise  from  free 
stream  velocity.  The  elastic  forces  are  internal  to  the  airframe.  Associated  with 
the  elastic  forces  are  inertial  forces  which  is  related  to  the  weight  distribution. 

Flutter  is  self-excited  vibrations  of  structures  caused  by  the  inability  of  the 
structure  to  dissipate  the  energy  received  from  the  air  stream.  Due  to  the  elasticity 
of  the  structure  and  the  external  force,  the  structure  gives  rise  to  torsional  and 
bending  motion.  As  the  structure  deflects,  a  new  geometry  is  presented  to  the  free 
air  stream,  thus  the  entire  cycle  of  bending  and  deflection  is  repeated.  This 
oscillating  motion,  according  to  linear  theory,  will  increase  exponentially  if  the 
structure  cannot  dissipate  the  external  energy. 

This  interaction  of  elastic,  aerodynamic  and  inertial  forces  on  structure  is 
usually  analyzed  by  taking  a  typical  section  of  the  wing  and  modeling  the 
torsional  and  bending  stifBiess  by  springs  (see  Figure  1).  At  a  specific  velocity, 
the  system  becomes  unstable. 
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Figure  1 .  Typical  section  in  flutter  analysis 
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The  two  equations  of  motion  that  describe  the  behavior  of  this  typical 
section  with  elastic  axis  at  the  mid-chord  are  derived  from  Hamilton's  principal 
and  are  shown  below 


Mg,-S,q,+K,q,  =  L  (2.1) 

-S,  q,  +  I.  q,  +K,q,=eL+  M.  (2.2) 

where 

M  mass/unit  length 

Sa  Static  moment/unit  length  (+  for  C.G.  aft  of  spring) 

Kh  Bending  stifftiess/unit  length 
K„  Torsional  stiffiiess/unit  length 

mass  moment  of  inertia/unit  length 
qi,q2  generalized  coordinates  for  bending  and  twisting  respectively 
L  Lift/unit  length 
Ma  Moment/unit  length 

2.2  Linear  Flutter  Analysis  -  Ouasi-steadv 

From  the  aerodynamic  theory,  lift  and  moment  can  be  formulated  in  the 
approximate  quasi-steady  manner  as 

Z,  =  7  p  c  a  (2.3) 

=  i  p  C/^  q,  (2.4) 

where  p  is  the  air  density  and  U  is  the  velocity.  The  angle  of  attack,  a  ,  can  be 
separated  into 

a  ^ 

a  =  or  <x  =  e-- 


(2.5) 
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Figure  2.  Elements  of  angle  of  attack 
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where  0  is  angle  of  attack  due  to  pitching  and  /i/U  is  angle  of  attack  due  to 
vertical  translation  of  the  aircraft  (see  Fig.  2).  Eqs.  (2.1)  and  (2.2)  then  become, 

L  =  qc  (q,  -  p  -  L’-d)  (2.6) 

eL  +  M,  =  eqcC,,  (g,  -p  +  qc^  C„  g,  +  M\t)  (2.7) 

where  is  the  dynamic  pressure,  p[/V2,  e  is  the  distance  between  elastic  axis  and 
the  mid  chord,  and  and  terms  describing  aerodynamic  forces  due  to 
disturbances.  Placing  Eqs.  (2.6)  and  (2.7)  into  (2.1)  gives 


M  q  + 


u 


(2.8) 


qSeC^ 


La 


U 


+  (K^-qSeCjq,  = 


(2.9) 


One  can  gain  much  insight  by  examining  the  system  without  damping. 
Setting  disturbance  loads  and  damping  terms  to  zero,  Eqs.  (2.8)  and  (2.9)  can  be 
written  in  a  matrix  form  as. 


Qa 

K-qSeC^ 


=  0 


(2.10) 


where  the  first  square  matrix  is  the  symmetric  mass  matrix  [mjj]  and  the  second 
square  matrix  is  an  unsymmetric  stiffiiess  matrix  [ky].  The  stiffiiess  matrix  varies 
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depending  on  the  dynamic  pressure.  It  is  this  unsymmetric  stiffiiess  matrix  which 
causes  system  instability. 

Assuming  sinusoidal  motion,  ,  Eq.  (2.10)  then  becomes, 


M  s  +  s  -¥ 


I  +  k 
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(2.11) 


For  stability,  examine  the  non-trivial  solution  by  setting  the  determinate  to 
zero,  gives. 


As*  +  Bs^  +C  =  0  (2.12) 

where  the  coefficients  are, 

A  =  MI^-Sl  (2.13) 

B  =  kuh  +  k,,M  +  k,,S^  (2.14) 

(2.15) 


The  roots  are  then  written  in  the  form  of  quadratic  formula, 


2  -B±Jb^-4AC 
=  - 2 - 

2  A 


(2.16) 


Stability  of  the  system  can  then  be  determined  by  examining  all  four  roots,  s, 
as  (I  increases.  At  low  values  of  q,  all  four  roots  are  imaginary  numbers.  At  a 
certain  velocity,  the  value  B^-4AC  becomes  negative,  as  a  result,  the  roots  become 
complex.  Dynamic  instability  occurs  at  the  point  where  the  real  part  of  any  root 
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becomes  positive.  In  addition,  at  a  specific  velocity,  two  of  the  roots  become 
zero.  This  signifies  static  instability. 

The  quasi-steady  method  of  describing  the  aerodynamic  forces  is  an 
extremely  crude  approximation,  but  it  gives  a  simple  physical  feel  for  the 
instability.  Note,  that  there  are  no  damping  terms  included  in  the  above  quasi¬ 
steady  approximation.  Due  to  the  inaccuracy,  analysis  of  events  even  with  quasi¬ 
steady  damping  may  result  in  an  incorrect  solution,  particularly  if  the  C.G. 
coincides  with  the  elastic  axis,  i.e.,  S„=0. 

2.3  Linear  Flutter  Analysis  -  Unsteady.  Fade  Appro^timants 

As  a  basis  for  comparison  for  the  full  nonlinear  flutter  analysis,  it  is 
beneficial  to  examine  the  linear,  small  amplitude  flutter  problem  with  zero  root 
angle  of  attack.  There  are  several  ways  to  solve  the  system  of  linear  flutter 
equations.  One  of  the  conventional  ways  is  Fade  Approximations.  Using  Eqs. 
(2.1)  and  (2.2),  the  unsteady  incompressible  aerodynamic  force  can  be  written  in 
the  form. 


Lea  =  Ti9b''[-h  +  U^-baQ 


+  2TipUbC{k) 


-^  +  t/e  +  6 


e 


(2.17) 


or  in  Laplace  Domain, 
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^EA  =  — T—  02% 


^  +  pQ  +  a  Q 
b 


+  ^^^47tC(/?) 


-p  +Q  +  \  a 


pe 


(2.18) 


where  b  is  the  semi-chord,  a  is  the  non-dimensional  parameter  defined  by  the 
distance  fi-om  the  mid-chord  to  the  elastic  axis  divided  by  b,  C(k)  is  the 
Theodorsen  function  and  p  is  pblU.  Using  the  single  lag  approximation  for  the 
Theodorsen  function  in  the  Lapace  Domain, 


C(p)  = 

p+.\5 

and  with  some  algebraic  manipulation,  Eq.(2.18)  becomes, 


(2.19) 


=  ^^b2% 


■p^-\Ap- 


.135^ 

^-h.15 


qU^ 

+  ^^b2% 


=2 


(\  V 

1-hl.l 

u  )\ 

-.9  +  . 135 


(\ 


p  +  2  +  - 


-a  \p 


P+.15 


e 


(2.20) 


further  simplification  yields. 
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Lsa  =• 


2 

pf/^ 


•®0/<  ■*■  _ 


^+.15 


^2bP  '*'  ^\bP  '*’  ^Ofl  _ 


^3bP 


P  +  A5 


h 

b 

e 


where, 


B,,  =  -2t: 

B,.  = 

B,,  =-2.2n 

B,.= 

;iu. 

II 

O 

^OB  ^ 

B,^  =  -.27071 

^35  ~ 

O  ' 

- a 


u 


-  a 


\ 

/ 


The  aerodynamic  moment  about  the  elastic  axis  (reference  line), 
found  in  a  similar  way, 


= 


p 

2 

p£/' 


^ic  P'  +  ^ir  P  +  ^oc 

^20  P'  +  ^1/2  P  + 


+  ^3C 


+  -^30 


P 

^+.15 

P 


A 

Z> 

0 


^+.15 


where, 


=  -2  a  7t 


^2D  ~ 


-2n 


1 

+  a 


V8 


J 


(2.21) 


(2.22) 


can  be 


(2.23) 
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=  -{2.2n)^^  +  a'' 


5,0  =  2  TT 


-  a 


■ 

-1  +  1.1 

J 

_ 

I2  )\ 

^oc=0 


^OD  =^n^  +  a  ^ 


(2.24) 


53,  =-(.270  71) 


O  ^ 

—  +  a 

U  j 


530  =  2  71  —  +  a 


■ 

11 

-.9+. 135 

J 

. 

I2  jj 

Here,  for  convenience  in  subsequent  nonlinear  analysis,  the  reference 
airforces  will  always  be  taken  at  the  quarter-chord,  so  that  a  =  -1/2.  This  gives 
the  following  simplified  form  of  the  coefficients, 


^2  A  =  -27: 
5„  =-2.271 
^0.  =0 

=--21n 


5.,  =0 
5o,  =  0 
53,  =0 


^2B  =  7: 

5,5  =  4.271 
Bob  =  471 
535  =  -1.5371 

5^0  =  -.7571 
5,0  =  -271 
Bqd  —  ^ 

B^D  — 


(2.25) 


Equations  (2.21)  and  (2.23)  would  then  express  L1/4  and  Mj/4  respectively, 
with  h  and  0  being  the  deflections  at  the  quarter-chord.  The  corresponding 
structural  coordinates  q,and  q2,  which  define  deflections  at  the  mid-chord  as 
shown  in  Fig.  1,  should  then  be  transferred  to  the  quarter-chord  using  the 
relations. 


I.  ^ 

0  =  ^3 


(2.26) 
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Inserting  Eqs.  (2.26)  into  the  unsteady,  incompressible  lift  equation  (2.21)  yields, 


I,  =  -pC/'  b 


+  B\aP  +  +  ^iA  - 


+  -pt/'Z) 
2 


^IB  + 


P+.15J  b 


(2.27) 


Inserting  Eqs.  (2.26)  into  the  moment  equation  (2.23)  similarly  yields. 


M,  = 


^IcP  ^\cP  ^OC  ^3i 


/?  +  .  15 


+  ^pt;v 


52£.  +r^:c  P  + 


^  ^\D  ■*■  2  ^ic  2  j 


+  1  Btn  +  T  ^3C  I  ^ 


2  ^^7p  +  .15j 


^2 


(2.28) 


Placing  these  into  the  basic  equations  of  motion  gives  the  right  hand  side  of  Eq. 
(2.1),  L,  as. 


I  =  -pt/'Z> 
2 

1 


^2aP  ^\aP  ^OA  ^2A  — 


p+.\5 \  b 
n  ,  D  P 


^2B  P  ^\B  P  ^OB  ^^B  — 


P+.15J 


(2.29) 


^2 


Similarly,  combining  Eqs.  (2.27)  and  (2.28)  to  form  the  right  hand  side  of 
Eq.  (2.2),  eL+M„,  gives, 
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|i  +  M.  =ipU^ 


1 


^2cP  ^\cP  ^OC  '*’  ^IC  — 


P+.15 


2 


^2dP  ^\dP  ■*■  ^OD  '*'  ^3D  ~ 


-5  P 


P  +  A5 


il 

b 

^2 


(2.30) 


where, 


Bu=B, 

Bib  =  B^b  +  -  5,^ 

4  =5„+l5, 

BjD  ~  B^^  +  —  B^  ^  ^/s  ^ 


(2.31) 


(/•  =  0,1,2,3) 


Then,  introducing  the  augmented  state  variables  Yg  defined  as. 


P 

p+.  15 
P 

p+.  15 


^2 


and,  reverting  back  to  time  domain,  yields. 


(2.32) 


B  -  -pU'^[B2^  qi  ■•r  +  B^^  q^  +  £3^  ] 

^  P  ^  ^[-^25  ^2  B^b  Q2  ■*■  Bqb  B^b  ^2  ] 


(2.33) 
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-L  + 
2 


—  ^pU  •^ic  ^1  ^OC  '*'  ^3C  ^  ] 

■*■  ^  P  ^  ^  [-^2/)  92  ■*■  ^ID  4 2  ■*■  ^OD  92  ^3D  ^2  ] 


(2.34) 


Y,+^{A5)Y,=q,  (2.35) 

b 

r,*^{.l5)Y,=g,  (2.36) 

b 

The  equations  of  motion  can  then  be  formed  by  placing  Eqs.  (2.33)  and 
(2.34)  into  the  right  hand  side  of  Eqs.  (2.1)  and  (2.2).  Together  with  the 
augmentated  state  equations  (2.35)  and  (2.36)  they  can  be  written  in  the  matrix 
form  as, 


'm* 

0 

o' 

f  •  ^ 

9 

■  0 

-M' 

0  ■ 

9 

0 

M' 

0 

< 

9 

>  + 

K' 

B' 

G’ 

•< 

9  ‘ 

0 

0 

I 

Y 

L  J 

0 

-I 

//• 

Y 

multiplying  by  inverse  of  first  matrix  gives, 


0 

0 


I 

-M'-'B' 

I 


0 

-H' 


9 

<  q  > 
Y 


(2.37) 


(2.38) 


where, 
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M*  =' 


M  -S^ 
-S^  / 


2  t/' 


5. 


2^ 


6  5- 
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-S. -ipA’B; 
-s. -ipi’s;,  /. -^pi'i 
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2D 
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-IpD^A 

2  U 
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5, 
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Z)5, 


15 


5,c  b^  5, 


ID 


--pUb^  5, 
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15 


ID 


K'  = 


0 

0  K 
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5, 
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bB, 


OB 


bB,,  b^  B. 


OD 


K, 


■  2 
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pU^B, 
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OA 


pU^B.^ 


3B 


--pWbB,, 
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K^--pW  b^B,, 


G* 


B. 
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bB. 
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-ipC/=i.B„ 


*B.,  ■ 
b"  B,r, . 

-|pC/’*B„ 

-jpC/' 


(2.39) 


(2.40) 


(2.41) 


(2.42) 
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.15 

0 


0 

.15 


(2.43) 


and  I  is  the  identity  matrix.  The  eigenvalues  of  the  matrix  in  Eq.  (2.38)  are 
solved  for  different  values  of  U.  In  this  study,  there  are  six  eigenvalues,  two  pairs 
of  complex  conjugates  and  two  real  roots.  These  roots  form  a  root  locus  plot. 
The  flutter  velocity  is  the  velocity  associated  with  the  complex  root  which 
changes  from  having  negative  real  value  to  positive  real  value.  In  addition,  the 
divergence  velocity  can  be  found  by  examining  the  two  real  roots.  One  of  the  real 
root  will  become  zero  at  the  divergence  velocity.  The  divergence  velocity  can 
also  be  found  from  the  matrix  K*.  The  velocity  which  makes  the  determinate  of 
K*  equal  to  zero  is  the  divergence  velocity.  The  results  using  the  values  listed  in 
Appendix  A  are  shown  in  Chapter  4. 


2.4  Linear  Flutter  Analysis  —  U-g  method 

Another  conventional  ways  to  analyze  linear  flutter  problem  is  the  U-g 
method.  As  expected,  this  method  starts  with  the  basic  equation  of  motion,  Eqs. 
(2.1)  and  (2.2).  The  forcing  terms  of  the  equations  are  formulated  by  assuming 
complex  sinusoidal  motion. 


e(0  =  ee'“' 


(2.44) 


The  aerodynamic  lift  for  2-dimensional,  incompressible  flow  is  then  given  by. 
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L  =  h  +  U mQ  +  (Si^  b a 0]e'“' 

+  2npUbC{k) 


i(iih  +  UQ  +  mb\  -  -  a  |e  |e'“' 


J_ 

2 


(2.45) 


where  C(k)  is  the  Theodorsen  function,  and  a  is  the  same  as  defined  in  Section 
2.3.  With  some  algebraic  manipulations,  Eq.  (2.40)  becomes. 


L 


(2.46) 


where, 


4  =  1- 


2iC{k) 


2C{k) 


L  —  fl  + 

0  7^2 


k^  ^  k 


\  +  2C{k) 


(\ 


- a 


(2.47) 

(2.48) 


Similarly,  the  aerodynamic  moment  can  be  found  by  assuming  complex 
sinusoidal  motion,  the  final  formulation  are  written. 


b* 


w, 


‘A  .  +  '"b  0 


(2.49) 


where. 


m. 


(2.50) 
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1  ^C(/:) 

-  +  a 
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—  a 

—  +  a 

\ 

2 

C{k)-\ 

V2  j 

) 

J 

(2.51) 


Again,  as  stated  in  Section  2.3,  the  reference  airforces  will  always  be  taken  at  the 
quarter-chord  for  later  nonlinear  convenience,  then,  a  =  -1/2.  The  corresponding 
strutural  coordinates  q,  and  q2,  which  define  deflections  at  the  mid-chord,  should 
also  be  transferred  to  the  quarter-chord  using  Eqs.  (2.27).  Eqs.  (2.41)  and  (2.44) 
then  becomes. 


L 


71  pO)^ 


1  ^ 

4+^4 

^2 

2  J 

(2.52) 


M  = 


1 


L 


2  j 


U  A) 


^2 


(2.53) 


Placing  a  =  -1/2  in  Eqs.  (2.42),  (2.43),  (2.45),  and  (2.46),  results  are  much 
simplified,  they  are. 


/  1 

L  =  1  -  z - 

*  k 

/  1  .  2C(^) 

®  2  e 
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+  f(l  +  2C(*)) 
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(2.54) 

(2.55) 

(2.56) 


/ 

J 


(2.57) 
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In  formulating  the  structural  terms,  instead  of  incorporating  viscous 
damping  terms,  structural  damping  terms  is  inserted  into  the  Eqs.  (2.1)  and  (2.2). 
The  resulting  equations  can  then  be  written  as. 


^  -  S,  ^2  +  K,{\-\-ig)  q^-  L  = 


-S,  +  h  ^2  +  (1  +  ig)  (I2  - 


(b 


-L+  M.  =0 


u 


(2.58) 

(2.59) 


Substituting  Eqs.  (2.41)  and  (2.44),  along  with  the  appropriate  sub-components 
and  Eqs.  (2.27),  into  Eqs.  (2.53)  and  (2.54),  yields  the  following  form  of  the 
equation  of  motion,  written  in  matrix  form. 


where 


(B„  -  K,  Z)  B„ 

B,,  {B,,  -  Z) 
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(2.60) 

(2.61) 

(2.62) 

(2.63) 

(2.64) 


(2.65) 
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The  solution  method  is  to  select  a  value  of  reduced  frequency,  and  solve  Eq. 
(2.31)  for  all  complex  eigenvalues  Zj.  For  each  value  of  Z,  the  corresponding 
structural  damping,  frequency  and  velocity  are  given  by. 


Im{Z}  1 

-  ;  CO  =  - 

Re{Z}  Tr^Iz} 


(2.66) 


The  procedure  is  repeated  for  several  values  of  the  reduced  frequency,  k, 
until  sufficient  number  of  data  points  have  been  generated  to  produce  a  smooth  U- 
g  diagram.  The  flutter  point  is  where  the  structural  damping,  g,  goes  to  zero.  The 
frequency  and  velocity  corresponding  to  that  point  is  the  flutter  frequency  and 
flutter  velocity  respectively. 

At  first  glance,  the  typical  section  analysis  seems  crude.  However,  the 
method  characterizes  the  entire  wing  fairly  accurate.  To  modify  the  typical 
section  for  an  actual  wing  with  varying  spanwise  deflection,  one  would  replace  q, 
and  q2  with  (l),(x)  q,  and  (|)2(x)  qj,  where  (t»,(x)  and  (t)2(x)  represent  the  first  bending 
and  first  torsion  modes  respectively.  Multiplying  Eq.  (2.53)  by  <j),(x),  Eq.  (2.54) 
by  <|>2(x)  and  integrating  over  the  span,  J  gives  (for  uniform  wing  properties), 

/,  M  q,  -  I,  q,  +  I,k,{\  +  ig)q,-\^,LcK  =  0  (2.67) 

-I,  q,  +  hI^q,^hkA\^ig)q.-\^2Kd5c  =  (i  (2.68) 

where  I„  I2,  and  I3  are  aerodynamic  integrals  defined  in  Appendix  A.  Dividing 
Eq.  (2.62)  by  I,  and  similarly,  divide  Eq.  (2.63)  by  I3  and  organize  in  matrix  form 
similar  to  Eq.  (2.55)  gives. 


35 


(B,,  -  KJ) 


The  determinate  is  then, 


(2.69) 


det[  ]^(B„-K,Z)(B,,-K,Z)-^B„B„  (2.70) 

The  accuracy  is  effected  by  the  ratio  I]  I  l^Iy  In  the  present  study  the  ratio  is 
0.92,  which  gives  only  a  small  correction  to  the  typical  section  results.  The 
results  of  U-g  analysis  is  shown  in  Chapter  4. 


CHAPTERS 


NONLINEAR  THEORY 


3.1  Nonlinearity  in  Aeroelasticity 

Nonlinearity  in  aeroelasticity  can  arise  from  either  the  structure  or 
aerodynamics.  Nonlinearities  related  to  structures  are  either  of  geometric  or 
material  origin.  Geometric  nonlinearities  depend  only  on  geometric  quantities, 
such  as  displacement  or  length,  such  as  in  large  deflection  of  beams  and  plates,  or 
large  displacement  of  helicopter  blades.  Material  nonlinearities  can  arise  from 
stiffiiess  properties  of  a  particular  structure.  On  the  other  hand,  nonlinearity  due 
to  aerodynamics  can  arise  from  stalling  and  the  attendant  flow  separation  of  an 
aircraft  wing. 

Nonlinear  limit  cycle  oscillations  can  originate  from  either  aerodynamic  or 
structural  nonlinearity  or  a  combination  of  both  [Ref  30].  A  simple  example  of  a 
generic  system  with  stiffening  springs  can  illustrate  the  concept  of  limit  cycle 
oscillation  (see  Fig.  3).  In  such  a  system,  if  assumption  of  linear  stiffiiess 
characteristic  is  prescribed,  then  the  oscillation  amplitude  would  grow 
exponentially.  However,  if  the  spring  behaves  nonlinearly,  as  the  oscillation 
penetrates  the  nonlinear  region,  the  effective  stiffiiess  may  increase  and  the 
oscillation  can  settle  down  to  a  constant  amplitude,  steady,  limit  cycle  oscillation. 
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The  present  investigation  will  concentrate  in  particular  on  the  aerodynamic 
nonlinearity.  The  including  of  structural  nonlinearities  is  discussed  briefly  in 
Appendix  F. 
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3.2  ONERA  Aerodynamic  Model 


The  aerodynamic  model  used  for  this  study  was  initially  developed  at  Office 
National  d' Etudes  et  de  Recherche  Aerospatiale  (ONERA)  by  Tran  &  Petot  [Ref. 
10]  and  by  Dat  &  Tran  [Ref  11].  This  aerodynamic  model  is  a  semi-empirical, 
unsteady,  nonlinear  model.  It  uses  quasi-linear,  small  amplitude  oscillation, 
experimental  data  to  predict  aerodynamic  forces  on  an  oscillating  airfoil  that 
experiences  dynamic  stall.  The  model  incorporates  the  Theodorsen  function  for 
linear  theory  by  inserting  a  single  lag  term  operating  on  the  linear  part  of  the 
airfoil's  static  force  curve.  The  model  also  includes  a  two  lag  term  on  the  stalling 
portion  of  the  airfoil's  static  force  curve. 

The  ONERA  model  was  later  investigated  by  Peters  [Ref  12]  to  distinguish 
between  the  angle  of  attack  due  to  pitching  and  angle  of  attack  due  to  vertical 
motion  of  the  aircraft.  The  final  form  of  the  ONERA  model  used  in  the  present 
study  are  shown  below. 


—  C^j  +  Cj2 


^Z\  ~  ^Z\  ^Z2  ^Zy 


*  ^ 

a  +  0 

y 


+  a^i^  I  a+  0 


•  «• 


Czy+  A,,  =  A., 

•*  • 

Cz2  +  r,  Cz2  +  rj  Qj  =  -r,  AC  L  -''3  ^  Cz |„ 


(3.1) 

(3.2) 

(3.3) 

(3.4) 


where  the  angle  of  attack  a,  is  separated  into  pitching  and  plunging  as  illustrated 
in  Fig.  2  and  the  non-dimensional  time  derivative  is  denoted  by, 


T  s 


b 


(3.5) 
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Figure  4.  Description  of  static  aerodynamic  curve 
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The  aerodynamic  force  curve  is  shown  in  Figure  4.  The  function  Cz  can 
represent  any  of  the  relevant  non-dimensional  aerodynamic  force  coefficients: 
Cl,  the  lift  coefficie'nts,  Cd,  the  drag  coefficient  or,  C^,  the  moment  coefficient. 
The  function  Cn  is  the  static  force  in  linear,  unstalled  flow,  extended  to  high 
angles  of  attack.  The  function  0^2  represents  the  nonlinear  aerodynamic  force 
contribution  to  the  total  Cz-  The  ACz  function  is  the  nonlinear  deviation  from  the 
extended  linear  force  curve.  In  addition,  ACz  also  the  input  which  determines 
the  magnitude  of  the  contribution  of  Cz2  to  the  total  aerodynamic  force.  The 
coefficients,  Szi,  Sz2,  §23,  A,],  X2,  rj,  r2,  r3  are  determined  empirically  and  are 
associated  with  the  appropriate  force  coefficient.  These  coefficients  are  listed  in 
Appendix  A, 

Equations  (3.2)  and  (3.3)  describes  the  linear  portion  of  the  model  where  Cl^ 
is  the  linear  circulatory  contribution  incorporating  aerodynamic  lag  due  to 
formation  of  tip  vortices.  The  function  ACz  is  defined  positive  for  a  decrease  in 
the  aerodynamic  force  beyond  stall  angle,  as  shown  in  Fig.  4.  The  general  form  of 
the  static  aerodynamic  force  curve  is  given  by. 


C,(a)  =  a^,a-AC,(a)  (3.6) 

where  a^z  is  the  slope  of  linear  aerodynamic  force.  In  general,  ACz  can  be 
approximated  in  any  manner  appropriate  for  the  particular  study.  In  the  present 
investigation,  the  objective  is  to  develop  a  simple,  nonlinear,  analytic  aeroelastic 
flutter  analysis.  The  function  ACz  is  therefore  described  by  a  single  straight  line 
fit  between  discrete  points  as  shown  in  Figure  5. 
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3.3  Harmonic  Balance  Method 

The  flutter  equations  and  related  components  have  all  been  stated  in 
differential  form.  It  remains  to  reduce  the  system  of  differential  equations  to 
algebraic  form  so  that  it  is  more  easily  solved  computationally.  By  assuming 
harmonic  motion,  the  harmonic  balance  method  is  a  powerful  way  to  transform  a 
system  of  nonlinear  differential  equations  to  a  simpler  algebraic  form. 

The  linear  and  nonlinear  aerodynamic  force  can  be  written  with  mean,  sine 
and  cosine  parts  as, 

Q  =  +  C^si  sin  kz  +  Qc,  cos  kz 

+  C^sj  sin2^x  +  C2C2  cos  2^1 

=  ^Zlo  ^ZIS  ^ZIC 

C22  ~  ^Z2o  ^Z2S\  sin  +  ^Z2C]  COS 

+  C22S2  sin2A-T  +  Q2C2  cos2A:t 


(3.7) 

(3.8) 

(3.9) 


3.3a  Linear  Aerodynamic  Forces 

As  a  preparation  of  nonlinear  flutter  analysis  and  the  harmonic  balance 
method,  harmonic  decomposition  of  the  aerodynamic  forces  were  necessary. 
First,  harmonic  motion  is  assumed  for  the  pitching  angle  and  deflection  at  the  1/4 
chord, 


where, 


sin^x  +  0^  cos^x 
sin^x  +  cos  kz 


(3.10) 

(3.11) 
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k  =  reduced  frequency  =  — 

T  =  non-dimensional  time  =  — 

b 

The  first  and  second  non-dimensional  time  derivative  of  Eqs.  (3.10)  and  (3.1 1)  are 
also  needed  for  subsequent  analysis,  they  are, 


Q  =  Q cos,  kz  -  Q s\n  kx  (3.12) 

•  « 

0  = sin  A:t  -  cos/:t  (3.13) 

h  =  h^k  cos  kx  -  h^k  sm  kx  (3.14) 

•* 

h  = -h^k^  sin  kx  -  h^k^  cos  kx  (3.15) 


The  angle  of  attack  is  a  combination  of  angle  of  attack  due  to  pitching  and 
angle  of  attack  due  to  plunging  as  stated  in  Eq.  (2.4).  Substituting  Eqs.  (3.10)  and 
(3.14)  into  Eq.  (2.4)  gives, 


a  =  +  a,  sin  kx  +  cos  kx 


(3.16) 


where 


a.  =  0. 


a^  =  e  +:ik 

S  ^  L 

D 


(3.17) 

(3.18) 

(3.19) 


Similarly,  Eq.  (3.3),  the  circulatory  part  of  the  linear  aerodynamics,  can  be  written 
in  the  form  of. 
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Qy  =  Qyo  +  Qy.  sin  kx  +  COS  kx  (3.20) 

Czy  =  k  COS  /:t  -  A:  sin  kx  (3.21) 

Substituting  Eqs.  (3.20)  and  (3.21)  into  the  left  hand  side  of  Eq.  (3.3),  and 
inserting  Eqs.  (3.10)  and  (3.11)  into  the  right  hand  side  of  the  same  equation. 
Matching  the  mean,  sine  and  cosine  terms  of  the  resulting  equation,  gives, 


c> 

CD 

II 

0* 

(3.22) 

=  F(k)  L,  -  G(k)  L, 

(3.23) 

=  G(i)  4  +  f(t)  4 

(3.24) 

where,  in  the  present  study,  the  F(k)  and  G(k)  functions  are  the  approximations  to 
the  Theodorsen  function,  C(k)  =  F(k)  +  iG(k),  namely. 


F{k)  = 


k%  + 

1]  +  e 


G{k)  = 


k  K  (K  -  0 
x\  +  e 


(3.25) 

(3.26) 


which  result  fi'om  using  the  single  lag  pole  approximation  of  the  generalized 
Theodorsen  function,  C{p)  =  {X2P  +  X^)l{p  +  X|),  with p  =  /k,  =  0.15,  A,2  =  0.55, 
and  where  other  intermediate  variables  are, 

4  =  ^02(^65  -  ^ y  -  ^ j  (3-27) 


(3.28) 
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Finally,  combining  the  harmonic  form  of  the  circulatory  terms,  Eqs.  (3.22), 
(3.23),  and  (3.24),  and  the  harmonic  form  of  the  remaining  linear  terms  in  Eq. 
(3.2)  gives  the  linear  portion  of  the  aerodynamic  coefficient  C^i  in  Eq.  (3.8)  as, 


C  =  C 

^Z\o  ^Zyo 


^Z]S  ~  ^Z\ 


^z\c  ^Z\ 


j 


+  C. 


ZyS 


A:  05  +  ^  -  ^^2  0c  +  5^3  ^  05  +  Q 

b  J 


(3.29) 

(3.30) 

(3.31) 


3.3b  Nonlinear  Aerodynamic  Forces 

For  the  nonlinear  portion  of  the  ONERA  aerodynamic  model,  similar 
harmonic  decomposition  procedure  was  performed.  As  stated  previously,  the 
aerodynamic  force  curve  is  modeled  by  a  single  break  point  approximation.  The 
nonlinear  deviation  from  the  extended  linear  force  curve,  AC^,  is  therefore 
parameterized  by  a  single  straight  line  with  the  form, 

AQ  =  (a  -  a, )  (3.32) 

where  ai  represents  the  stall  angle.  Equation  (3.32)  is  only  valid  in  the  stalled 
region  described  in  Fig.  5.  The  manipulations  of  the  equations,  can  be  further 
simplified  if  the  angle  of  attack  is  put  in  the  form  where  it  is  purely  sinusoidal, 


a  =  a„  +  a^,  sin  ({) 


(3.33) 


47 


where, 


a,.  =  yjal  +  al 
<j)  =  ^ 


(3.34) 

(3.35) 


C  •  -I 

4  =  Sin  — 

a.. 


(3.36) 


The  formula  for  AC^  can  also  be  written  in  the  form  of  sinusoidal  components, 


AQ  =  +  AQj,  sin(j>  +  AC2C2  cos2(f> 


(3.37) 


Note  that,  in  the  above  formulation,  there  is  only  sine  term  associated  with  the 
first  harmonic.  This  is  due  to  the  fact  that  AC^  is  a  single-valued  function  of  a. 
The  two  functions  are  always  in  phase  with  each  other,  therefore  there  is  no 
cosine  term  in  Eq.  (3.37).  Similarly,  there  is  only  cosine  terms  associated  with  the 
second  harmonic.  Substituting  Eq.  (3.33)  into  Eq.  (3.32),  Fourier  analysis  then 
gives  the  relations  for  AC^o,  AC^si,  and  AC2C2. 


=  7  (a„  +  Ot^.  sin  ij)  -  a, )  t/cf) 

(3.38) 

=  7  A  (<^c  +  sin  <t)  -  a, )  sin  (j)  d<l> 

(3.39) 

^Qc2  =  7  *1  (cto  +  ciy  sin  <t>  -  a, )  cos  2<|>  4 

(3.40) 

After  some  algebraic  manipulation,  the  equations  become. 


AC 


Zo 


Z>1  ay 
n 


+  cos  ((), 


(3.41) 
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= 


_  b.CLy 


% 


<l>i 


71 


^sin2(t>, 


(3.42) 


AC 

71 


cos(|),  -  ^  sin3(|), 


(3.43) 


The  parameter  (j)i  a  non-dimensional  time  variable  associated  with  angle  of  attack, 
has  relationships  shown  in  Figure  6  and  by  the  equations  below, 


(j),  =  sin  '  5  if  -  1  <  5  <  1 


(3.44) 


<1>,  = 


f 

■|  if  5<-l 


(3.45) 


where, 


6  =  ^^ 


a.. 


(3.46) 


For  those  regions  where  (}>,  is  undetermined,  it  takes  on  the  values  of  ±  Till,  as 
described  in  Eq.  (3.46).  These  values  are  arbitrarily  set  so  that  the  limits  of 
integration  are  correct  in  the  Fourier  analysis  in  Eqs.  (3.38),  (3.39)  and  (3.40).  In 
the  present  analysis,  the  flutter  oscillation  amplitude  may  reach  the  negative  side 
of  the  aerodynamic  force  curve,  see  Appendix  D  for  the  formulation  of  the 
symmetric  force  curve. 


With  all  the  components  needed  to  formulate  the  nonlinear  aerodynamic 
force  in  harmonic  form,  it  is  just  a  matter  of  algebra.  Substituting  Eq.  (3.9)  and  its 
non-dimensional  time  derivatives  into  left  hand  side  of  Eq.  (3.4).  Similarly, 
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inserting  Eq.  (3.37)  into  the  right  hand  side  of  the  Eq.  (3.4)  and  express  this  right 
hand  side  as, 


R.  H .  S  =  sin  <|)  +  cos (j) 

+  sin  2(j)  +  cos2<t) 

where 


R,=-r2  AQ,  (3.48) 

R,,  =  -r,  AC„,  (3.49) 

Rc\  =  “'3  ^  AC31  (3.50) 

R,,=2r,k^C,,,  (3.51) 

=  -r,  AC^,  (3.52) 


The  harmonic  components  of  nonlinear  portion  of  aerodynamic  force,  €22?  aro 
found  by  matching  terms  of  both  sides  of  Eq.  (3.4)  then  gives. 


where. 


C22  —  Bxo  ^ZS\  ^  ^ZC\  cos  (j) 

+  5^52  sin24»  +  B2C2  cos 24) 


B 

Bzo  ~ 


p  —  •^i  ^s\  ^2  ^1 

+  Kl 


^ZCl  ~ 


—  ^1  ^1  ^2  ^s\ 

Kf  +  Kl 


Bzs2  — 


ATj  R, 


C2 


Kl  +  Kl 


(3.53) 


(3.54) 

(3.55) 

(3.56) 


(3.57) 
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D  _  ^3  ^C2  ^4  ^S2 

Kl  +  K] 


where  the  intermediate  variables  are, 


(3.58) 


II 

1 

KJ 

(3.59) 

K,=r,k 

(3.60) 

K2=r2-Ae 

(3.61) 

K,  =2r,k 

(3.62) 

The  above  results  are  formulated  with  the  angle  of  attack  in  the  purely 
sinusoidal  form,  where  all  trigonometric  functions  are  taken  with  respect  to  the 
non-dimensional  parameter  (().  To  make  subsequent  analysis  simpler  and 
consistent  with  the  linear  portion  of  the  aerodynamic  force,  inserting  ^  ^ 

from  Eq.  (3.35)  into  Eqs.  (3.33)  and  (3.53),  expanding  out,  gives, 

a  =  a„  +  a^.  cos  ^  sin  kx  +  sin  ^  cos  kx  (3.63) 


C22  =  B2^  +  B2SI  cos  4  sin  sin  ^  cos  kx 

B2CI  cos  ^  cos  kx  -  B2C1  sin  ^  sin  kx 
Bzsi  cos2^  sin2A:T  +  B2S2  sin2^  cos2A:x 
B2C2  cos 2 4  coslkx  -  B2C2  sin 2 4  sin2^T 


(3.64) 


Introducing, 


e  “s 

cos  q  =  — 


(3.65) 


52 


then,  Eqs.  (3.63)  and  (3.64)  become, 


a  =  a„  +  sin  kx  +  cos  kx 

—  C220  ^Z2si  sinA:T  +  C^jci  cos^t 

+  Q252  sin2/tT  +  C22C2  coslkx 


where. 


^Z2o  ~  ^Zo 


n  —  R  —  R 

'^Z2S\  ~  ^ZS\  _  ^ZC\ 


a. 


a. 


r  -  R  ft 

'-'Z2C1  “  "ZSl  "ZCl 

U,y 


^Z2S2  ~  ^ZS2 


'^S 

yal 


alj 


-  B. 


ZC2 


2^^ 
OLy  CLy 


C22C2  =  B2S2l^^  ^  B. 


ay  ay 


ZC2 


yOj, 


2 


a 


V  J 


(3.66) 

(3.67) 

(3.68) 

(3.69) 

(3.70) 

(3.71) 

(3.72) 

(3.73) 


53 


3.3c  Summary  of  Complete  Aerodynamic  Forces 

The  formulation  of  the  ONERA  aerodynamic  forces  with  harmonic  balance 
method  has  been  completed.  This  section  summarizes  the  results  in  preparation 
for  the  assembly  of  full  flutter  equations.  The  complete  aerodynamic  force 
coefficient  Cz  is  written  as, 


Q  =  Qo  +  Qsi  sin  kx  +  cos  kx 

+  C2S2  sin2A:x  +  cos2^t 

where  each  of  the  harmonic  components  are, 

C  =  C  +  C 

'-'Zo  Wlo  ^  '^Z2o 

c  =  c  +  c 

^ZSl  WlS  ^  ^251 

c  =  c  +  c 

^ZCl  ^Z]C  ^  '“22CI 

c  =  c 

^ZS2  '^Z2S2 

c  =  c 

^ZC2  ^Z2C2 

linear  aerodynamic  coefficient,  Czi,  was  found  to  be: 


^zi  ^zio  ^zis  sin  cos 


where. 


C  =  C 

Wlo  Wyo 


^ZIS  ~  ■^Zl 


^ZIC  “  ^Zl 


h  \ 

-ke,+^k^ 

{  '  b  ) 


h  \ 

'  b  J 


^Z2  ^  ^Z3  ^  ®C  ^ZyS 


^Z2  ®C  ^  ^Z3  ^  ^ZyC 
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and  the  intermediate  variables  are, 


^Zyo  ~  ^oZ 


=  F(k)  L,  -  G(k)  L, 
=  G(k)  L,  +  F(k)  L, 


F(k)  = 


Gik)  = 


X]  +  e 

kx,  (?i,-l) 

+  k^ 


k's  ~  ^oZ 


k'c  ~  ^oZ 


he 


e-k-^-kQr 

S  ^  C 


h  ^ 

+  k  —  +  k  Q. 

^  b  'j 


The  nonlinear  aerodynamic  coefficient  is, 


C22  —  G220  ^Z2S\  k"^  +  C 


Z2C1 


^Z2S2  2^T  +  C22, 


=  B^ 

~  -^ZSI 

Bzc\ 

ac 

ay 

ay 

“  Bzs\ 

“c 

+  Bzc\ 

ay 

ay 

coskx 
,,  coslkx 


where, 
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^Z2S2  ~  ^ZS2 


21 

yal 


2  \ 

£ 

2 

y  ) 


-B  2^^ 

^ZC2  ^ 

a^.  a^. 


C  —El  2L  2:1.  +  f) 

^Z2C2  ~~  ^  "r 


'ZS2  ‘ 


CLy  ay 


ZC2 


2 


a 


2  \ 


ya;  a].j 


and, 


^zs\  ”■ 


^ZCl  ^ 


^ZS2  “ 


^ZC2  “ 


^S\ 

+ 

^2 

^C\ 

k; 

+ 

Kl 

K,Rc, 

^2 

^S\ 

K! 

+ 

^3 

^S2 

+ 

K, 

^C2 

K] 

+  ^4 

^3 

Bc2 

— 

^S2 

Kl 

+ 

K] 

R,  =  AC,^ 


"ZSl 


/?5,  =  —f’2  AC^; 

Rc\  ~  ~^3  ^  ^^ZSl 

Rg2  —  2  r^k  AC2(-2 

Rc2  =  ~^2  ^^ZC2 


K,=r2-  e 
K2=r,k 
=  r2  -  4 
K,=2r,k 
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n 


«,  -g. 

a,. 


n 


V 

^<t>,  --J  +  cosf 


AC 

71 


^  1  . 


<j),  -- 
2  J 


--sin(2(f),) 


^'^ZC2  _ 

7t 


-  ^  cos<t),  -  ^  sin  3(t>, 


All  components  of  the  aerodynamic  forces  has  been  transformed  from 
differential  form  to  algebraic  form  using  the  harmonic  balance  method.  In  this 
form,  numerical  solution  method  can  be  easily  implemented.  Before  proceeding 
to  assemble  the  full  3-dimensional  flutter  equation  several  more  steps  remain: 
First,  transform  the  general  flutter  equations  to  full  3-dimensional  flutter  equation. 
Second,  transfer  the  aerodynamic  force  formulations  to  structural  coordinate. 
Finally,  transform  the  structural  terms  of  the  resulting  flutter  equations  to 
algebraic  form  and  assemble  the  full  flutter  equations. 


3.4  Nonlinear  Flutter  Analysis 


3.4a  Formulation  of  3-Dimensional  Flutter  Equations 

A  general  Rayleigh-Ritz  formulation  is  used  to  approximate  plate  deflection 
and  perform  modal  analysis.  The  Rayleigh-Ritz  analysis  begins  by  assuming  a 
deflection  shape  for  the  structure.  If  only  out-of-plane  deflection,  w,  and 
rotational  displacement,  a  are  allowed,  the  general  deflections,  written  with 
generalized  coordinates  are. 
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^=tySx,y)q^  (3.77) 

/  =  ! 

a  =  e«  +  I§^^,  (3.78) 

where  y,  (x,>^)  represent  the  non-dimensional  mode  shapes  of  the  i-th  mode;  qj 
denote  the  generalized  displacements,  or  modal  amplitudes,  of  the  i-th  mode;  n  is 
the  number  of  mode  shapes;  and  0r  is  the  root  angle  of  attack. 

For  simplicity,  it  is  further  assumed  that  the  mode  shapes  are  separable  in  the 
chord-wise  and  span-wise  directions.  The  mode  shapes  can  then  be  written  in  the 
form, 

Y,  (x,y)  =  (|).^(x);  for  bending  modes  (3.79) 

y  . 

y,(x,>')  =  —  <j)^„(x);  for  torsion  modes  (3.80) 

c 

Note,  in  the  present  investigation,  to  simplify  mathematical  work  and  still 
sufficiently  describe  the  deflection  of  the  wing  in  flutter  analysis,  only  beam  out- 
of-plane  bending  modes  and  beam  torsion  modes  were  represented.  Previous 
studies  have  used  varying  degree  of  complexity  in  terms  of  selecting  mode 
shapes.  Landsberger  and  Dugundji  [Ref.  25]  used  simplified  sinusoidal  torsional 
mode  shapes,  and  the  works  by  Dunn  [Refs.  13  and  14]  incorporated  mode  shapes 
with  higher  degree  of  complexity.  It  was  found,  however,  that  using  of  mode 
shapes  which  did  not  meet  the  cantilevered  root  warping  condition,  somewhat 
affected  the  Rayleigh-Ritz  prediction  of  modal  deflection.  Therefore,  more 
complex  torsional  modes,  which  incorporated  the  root  warping  effect  should  be 
considered  if  highly  accurate  description  of  wing  shape  is  required  [Refs.  26  and 
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27].  In  present  analysis,  root  warping  terms  will  not  be  incorporated  into  the 
mode  shape.  Furthermore,  since  the  final  stall  flutter  problem  is  expected  to  yield 
a  single  degree  of  fi-eedom  motion  in  either  the  first  torsional  or  first  bending 
mode,  it  was  decided  to  only  model  the  first  torsion  and  first  bending  modes.  The 
selected  mode  shapes  are  listed  below.  The  out-of-plane  bending  mode  is  written 
in  the  form. 


<t>/,(^)  =  cosh 


e, - 


1) 


cos 


X 

^'7, 


a, 

sinh 

(  A 

-  sin 

r 

-11 

1  '  i) 

1 

e,  =  p,  ;:  =  1.87510 


sinhe,  -  sine, 

a,  = - ! - L  =  0.72664 

COShE,  +  COSE, 


p,  =  0. 59686 
The  torsion  mode  is. 


<j>„(x)  =  sin 


y2l  j 


(3.81) 


(3.82) 

(3.83) 

(3.84) 


(3.85) 


Applying  the  Rayleigh-Ritz  energy  method  to  find  the  entries  in  the  mass 
and  stiffiiess  matrix,  where  the  kinetic  energy,  T,  for  the  system  is. 
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T  =  mw^  dydx 

1  i  J 

(3.86) 

M^j  =  JJwy.y^  dydx 

(3.87) 

The  internal  energy  formulation  gives. 

^  .mg  +  2£>,2W„w^+---]  dxdy 

=  ^11^.9,  K, 

Z  /  J 

(3.88) 

Placing  the  kinetic  and  potential  energy  terms  into  Lagrange's  equations,  results  in 

a  set  of  ordinary  differential  equations  of  motion, 

[M]{q}  *[K]{q)  =  [Q] 

(3.89) 

where  [M]  and  [K]  are  mass  and  stiffiiess  matrices  respectively. 

{Q}  is  derived  from  the  expression  of  the  virtual  work , 

The  modal  force 

hW  =  dy dx 

=  a 

1 

(3.90) 

Q.  =\ly,hpdydx 

(3.91) 

From  the  above  formulation,  the  entries  of  the  mass  matrix  then  are, 
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=  MU, 

(3.92) 

M„  =  /, 

(3.93) 

0 

II 

II 

(N 

(3.94) 

where  I,  =  \  dx,  \  dx ,  rg  denotes  the  radius  of  gyration,  /  represents 

the  length  of  the  wing,  and  c  is  the  chord  length. 

Again,  here  M  is  mass  per  unit 

length.  The  entries  in  the  stiffhess  matrix  [K]  are. 

(3.95) 

K  -  ^  Ao  j 

(3.96) 

where  D,  i  and  denote  the  flexural  rigidity  modulus  data  for  the  composite 
wing.  The  modal  forces  are  given  by  Eqs.  (3.90)  and  (3.91), 


Qi  =  Jo  dx 

Qi  =  —  j'  M  dx 
«  ^  Jo  aX 


(3.97) 

(3.98) 


where  L1/2  and  are  lift  and  moment  at  the  mid-chord.  These  aerodynamic 
forces  are  related  to  the  previous  formulation,  lift  and  moment  at  quarter-chord, 
as. 


~  h 


(3.99) 
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c  . 


(3.100) 


Now,  to  make  Eq.  (3.89)  easier  for  subsequent  analysis,  non- 
dimensionalization  of  the  equations  are  needed.  This  could  be  accomplished  by 
introducing  and  t,  where  t  represents  the  non-dimensional  time  as  stated  in  Eq. 
(3.5),  and  q-  denotes  the  non-dimensionalized,  generalized  coordinates  where. 


(3.101) 


and  from  Eq.  (3.5)  the  normal  time  derivative  can  then  be  written  in  the  form. 


—  =  -  IL(\ 

dt  b  dx  b 


(3.102) 


Substituting  Eqs.  (3.101)  and  (3.102)  into  Eq.  (3.89).  With  some  algebraic 
manipulation,  Eq.  (3.89)  can  then  be  written  as. 


PTt/,  q^+CTklq, 


~  Jo  *1^/.  ^ 

=  Jo  ^ 


(3.103) 

(3.104) 


where,  n  =  /  o)^,  the  ratio  of  the  bending  to  torsional  frequency,  k„  =  (D„b/U  is 

the  torsional  reduced  frequency,  and  the  rest  of  the  parameters  are. 


M 

npb^ 


(3.105) 
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(3.106) 


and  Ij  and  are  aerodynamic  integrals  for  the  selected  mode  shape  (see 
Appendix  A). 


3.4b  Coordinate  Transfer  of  Aerodynamic  Forces 

At  present  moment,  all  aerodynamic  forces  has  been  derived  about  the 
quarter-chord  as  a  matter  of  convention.  However,  the  structural  components  in 
this  study  have  been  formulated  with  the  longitudinal  axis  placed  at  the  mid¬ 
chord,  again  as  a  matter  of  convention  in  structural  analysis.  The  objective  of 
present  study  is  to  investigate  the  stalled  flutter  behavior  of  aircraft  structures.  It 
was  therefore  decided  that,  the  aerodynamic  force  parameters  should  be 
transferred  to  conform  to  structure  coordinates. 

Transferring  of  the  aerodynamic  forces  starts  by  converting  the  basic 
parameters,  pitching  angle  and  vertical  translation.  In  structural  coordinates,  the 
pitching  angle  and  vertical  translation  are, 

e,  +  (3.107) 

2  2 

(3.108) 

2 

The  pitching  angle  does  not  change  when  transferring  coordinates.  However,  in 
order  to  express  the  vertical  translation  about  the  quarter-chord  in  structural 
coordinates,  additional  terms  need  to  be  added. 
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A.  =*1+^6,  (3.109) 

4  2/2 

substituting  Eqs.  (3.108)  and  (3.109)  into  Eq.  (3.107)  and  denoting  Eq.  (3.108)  in 
the  quarter-chord  subscript,  then  gives, 


J 


0,  =  0 


+|<l>«^2 


(3.110) 

(3.111) 


Again,  introducing  harmonic  motion,  now  in  terms  of  the  non-dimensionalized  q, 
Qi  =  lo  +  ^is  sin  kx  +  cos  kx  (3.112) 


Then,  substituting  Eq.  (3.112)  into  Eqs.  (3.110)  and  (3.111),  result  in  the 
harmonic  form  of  vertical  motion  and  pitching  angle, 

K  -  K  sin  kx  +  h(.  cos  kx  (3.113) 

4 

where. 


A.  =  i( 


A  -  ix  - 
y  +  <!>*  ^10  +4‘l>a  ^20 


A  -  lx  -  1 

TA  Q\S  .  9o  ^IS 

4  y' 


X  -  1x0 

TA  ^ic  ^  Yo  ^2C 

4  y 


(3.114) 

(3.115) 

(3.116) 


and. 
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0^  =  0^  +  0^  sin /:t  +  cos  A:t  (3.117) 

4 

where, 

0.  =  0« +  (3.118) 

=  (3.119) 

With  all  the  building  blocks  at  hand,  transferring  of  the  aerodynamic  force 
coefficients  can  now  proceed.  Obtain  the  linear  aerodynamic  force  coefficient, 
Cz\,  by  inserting  the  harmonic  components  of  vertical  displacement,  Eqs.  (3.1 14), 
(3.115)  and  (3.116),  along  with  the  harmonic  components  of  the  pitching  angle, 
Eqs.  (3.118),  (3.119)  and  (3.120),  into  Eqs.  (3.29),  (3.30)  and  (3.31),  and  their 
appropriate  sub-components,  gives, 


"Zlo 


=  a 


oZ 


y 


^Z\S  ~  ^hi^Z\  ^\S  ^Z2  7ir)  **"  *l^a(^Z3^25  ^ZA  ^2C  ) 

^ZIC  ~  4**(~'^Z2  ^\S  ^Z1  ^1(  )  4*a(~^Z4  ^2S  ^Z3  ^2C  ) 


(3.121) 

(3.122) 

(3.123) 


The  above  formulations  seems  complicated.  However,  there  are  only  four 
components  involved  in  the  entire  formulation,  namely,  Eij,  E22,  ^zi»  ^Z4- 
Notice  the  symmetries  that  exist  between  the  Czjg  and  C^jc-  Each  of  the  sine 
components  of  C21S  is  that  of  the  cosine  components  in  Czicj  each  of  the 
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cosine  components  of  Cz,s  is  the  negative  value  of  the  sine  components  in  Czi^. 


These  components  are  written  as, 

+a^^G{k)k  (3.124) 

E,,=a„,F{k)k  (3.125) 

Ezi  —  ~  ^  (3.126) 

Ez4  ~  ~  2^21  ^  ~ ^  ~ '^^oz  ^  E(k)  — —0^2  G{k)  (3.127) 


The  nonlinear  aerodynamic  force  coefficient,  C22  can  be  transferred  to  the 
structural  coordinate  in  the  similar  fashion.  The  angle  of  attack,  a,  in  the 
structural  coordinate  is  needed  in  the  transferring  of  nonlinear  force  coefficient. 
By  inserting  Eqs.  (3.113)  and  (3.1 17)  and  their  non-dimensional  time  derivatives 
into  Eqs.  (3.17),  (3.18),  and  (3.19),  the  angle  of  attack  in  structural  coordinates  is 
then, 

sin^x  +  a,-  cos^x  (3.128) 

where, 

(3.129) 

“s  =  I  <l>a  ^25  +  <1>*  ^  ^IC  +  ^  K  ^^2C  (3.130) 

-Kkq^s (3.131) 

The  above  formulations  can  then  be  placed  into  Eq.  (3.68)  and  its  appropriate  sub¬ 
components  to  form  the  transferred  €22-  The  full  expression  will  not  be  listed 
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here  due  to  its  complexity.  When  solving  for  the  full  flutter  solutions,  is 
programmed  in  the  computer  code  as  it  appears  in  Eq.  (3.68).  The  code,  hoAvever, 
is  specified  to  place  the  transferred  values  of  a,  0  and  h  into  the  appropriate 
components  as  required. 

3.4c  Assembling  of  Flutter  Equations 

Once  the  aerodynamic  forces  are  set  up  for  this  analysis,  it  still  remains  to 
assemble  the  full  flutter  equations.  Employing  the  harmonic  balance  method, 
forcing  terms  of  the  full  flutter  Eqs.,  (3.103)  and  (3.104),  can  be  expressed  in  their 
harmonic  elements.  Recall  Eqs.  (3.99)  and  (3.100),  where  lift  and  moment 
coefficients  in  structural  coordinates  are  related  to  their  counter  part  in 
aerodynamic  coordinates.  Inserting  those  formulations  with  their  appropriate 
elements,  namely  Eqs.  (3.121),  (3.122)  and  (3.123)  and  Eqs.  (3.68)  with 
transferred  values  of  a,  into  the  forcing  terms,  then  gives, 

Jo  4»>.  ^  =Jd  +  Q.5.  sin  kx  +  coskx]  dx 

yjAjZ) 

+  Jd  [Cl2o  +  Q251  sin  kx  +  Qjc,  cos  kx]  dx 

Jo  ‘t>a  ^  =Jd  *t'a  \^M\o  +  ^  QlO  j  +  ^0/151  +  Ql5l  j  ^  QlCI  j  ^OS  A:T  dx 

+  /d  <1*0  ^  Cu2o  +  ^  Qzo  j  ^  Cuis\  +  ^  ^L2s\  j  sin  Art  +  ^  C*^2ci  +  ^  Q2C1  ^  cos  kx  dx 

(3.133) 

Notice  that  Eqs.  (3.132)  and  (3.133)  contains  only  first  harmonic  terms.  In 
the  present  study,  it  is  not  the  intention  to  model  the  fine  details  of  the  nonlinear 
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motion.  Therefore,  harmonic  balance  method  was  employed,  which  is  suitable  for 
describing  the  gross  aspects  of  the  nonlinear  problem.  This  idea  is  further 
amplified  by  the  use  of  a  semi-empirical  aerodynamic  model,  which  ignores  the 
fine  details  of  the  fluid  flow.  It  was  therefore  decided  to  use  only  first  harmonic 
terms  for  the  total  flutter  analysis,  while  the  aerodynamic  analysis  by  itself, 
namely,  fitting  of  the  aerodynamic  model  (see  Appendix  B),  can  use  several 
harmonics.  Furthermore,  the  angle  of  attack  at  75%  span  is  taken  to  be  the 
average  angle  of  attack  for  the  entire  wing  span.  This  assumption  greatly 
simplifies  the  formulation  and  is  implemented  in  the  nonlinear  portion  of  the 
aerodynamic  forces  by  evaluating  all  mode  shapes  at  the  specified  span  location. 

With  additional  algebra  work,  Eqs.  (3.132)  and  (3.133)  expand  to  the  form. 


=  A  0,  +  A  ^9.. 

[-^1  (^il  ^L2  ^\C  )  '*'  A  (^i3  ^25  ■^i4  Q2C  )] 

+  [^l(-^i2  +  ^ic)  +  A(-^i4  ^2S  +  -£^i3^2c)]cOSA:T 


— 

“ 

+  ^4  + 

h 

D  n 

"i51  ^LC\ 

y  CLy  ay  ) 

^LS\  ^LC\ 


L  V 


a 


a 


V  J 


cos 


(3.134) 
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Jo  ^  ^  0.  +  A  -f^2o 

[-^2  i^T\  ^\S  ^  ^T2  9\C  )  "*"  -^3  (^r3  ^2S  ^TA  ^2C  )] 

+  [^2  (”  ^T2  ^\S  ^T\  ^\C  )  ^2S  ^T3  ^2C  )] 


+  /j  + 


” 

f 

Is 

^751 

n 

•°7C1 

O-y  J 

sin  kz 


+ 


■^751  "''  ^TC\ 


\ 


a. 


a 


y 


cos  A:t 


(3.135) 


where  variables  with  a  subscript  L  denote  components  associated  with  lift 
coefficient  and  I's  are  the  aerodynamic  integrals.  Components  with  a  subscript  T 
represent  terms  associated  with  transferred  values  of  moment  coefficients.  They 
are  written  as. 


^oT  ~  ^oM  ^ 

E 

E  -  E  +  ^ 

•^ri  ^M\  ^  ^ 

E 

E  -  E  + 

^T2  ^M2  ^  ^ 

E 

E  -  E  + 

^M3  ^  ^ 

E 

Er,=E^,+^ 

Bu=B^+^ 


(3.136) 

(3.137) 

(3.138) 

(3.139) 

(3.140) 

(3.141) 


4 


(3.142) 
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^Tc\  ~  ^Mc\  ^  (3.143) 

Finally  assembling  the  foil  flutter  equations  by  introducing  Eq.  (3.112)  into 
left  hand  sides  of  Eqs.  (3.103)  and  (3.104),  the  differential  form  of  the  equation  of 
motion.  Similarly,  substituting  Eqs.  (3.134)  and  (3.135)  into  the  right  hand  sides 
of  the  equations  of  motion.  Matching  terms  of  both  sides  would  give  a  system  of 
six  coupled,  nonlinear  equations.  They  are  of  the  form, 


p  7r  /,  Q  =  I,  0^  +  ^  5 


Lo 


h  K  ^20  =  ^  ^^0 


(3.144) 

(3.145) 


p  71 7,  (q  k  )  —  /|  4'i5  ^L2  ^\c ) 

A  i^L3  ^IS  ■*'  ^2C  )  A 


^LS]  ^LC\ 


a 


a,  j 


p  71 7j  (q  k^  k  )  q^^  ~  (  ^12  ^\s  ^ L\  ^ic ) 

/ 

A  (“-^14  ^2S  ^L2  ^2C  )  "^  -^4  ^LC\ 


a. 


(3.146) 


R  ^ 

^LS\ 

ay  J 


(3.147) 


^  -^3  (^a  ^  )^25  “  -^2  (-^ri  ^'iS  ^T2  ^\c) 


+  /j  (£;-3  ^2S  •^7'4  92c)  :^5 


D  _  D  ^ 

•“TSI  "TCI 

a,  a 


V 


r  y 


(3.148) 
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J'-a  h  (K  -  (-£n  +  E„  q^,) 


"^3  (  •^r4  ^2S  ^2c)  ■*■ 


B  ~  +  B  ^ 

■°7C1  ..  +  ^TSi  - 


a. 


a. 


(3.149) 

Equations  (3.144)  through  (3.149)  are  the  full  flutter  equations  which  this 
study  is  based  on.  Notice,  again,  the  symmetries  that  exist  between  Eqs.  (3.146) 
and  (3.147)  .  In  both  the  linear  and  nonlinear  portion  of  the  aerodynamic  force 
terms,  the  coefficient  appearing  with  the  sine  terms  in  Eq.(3.146)  are  same  as  the 
cosine  term  coefficients  in  Eq.  (3.147).  The  cosine  term  coefficients  are  the 
negative  of  the  coefficients  appearing  with  sine  terms  in  Eq.  (3.147).  Similar 
symmetries  also  exist  between  Eqs.  (3.148)  and  (3.149).  These  symmetries  reflect 
the  fact  that,  if  complex  sinusoidal  motion  is  assumed  for  the  generalized 
coordinates,  these  six  equations  could  be  combined  and  form  a  system  of  three 
complex  equations.  This  system  of  nonlinear  Rayleigh-Ritz  aeroelastic  equations 
can  then  be  solved  by  using  Newton-Raphson  iteration  method  (see  Appendix  E). 

The  theoiy  described  in  all  previous  sections  was  implemented  using 
FORTRAN  code  with  Microsoft  FORTRAN™  compiler  Ver.  5.10  on  a  IBM 
compatible  386,  25  MHz  personal  computer.  The  source  code  of  all  programs 
written  are  listed  in  Appendix  F.  The  Newton-Raphson  solver  employs  a  finite 
difference  scheme  for  all  numerical  differentiation.  The  libraiy  routines  of 

LUDCMP  and  LUBKSB  [Ref  40],  were  also  employed  “to  perform  matrix 
inversion. 


CHAPTER  4 


RESULTS 


A  summary  of  linear  flutter  analysis  results  are  shown  in  Table  1.  The 
results  are  also  shown  in  graphics  form  in  Figs.  7  through  9.  All  linear  analysis 
results  were  generated  using  the  physical  properties  listed  in  Appendix  A,  Table  4. 
The  analytic  results  for  linear  flutter  thoery  using  the  Fade  Approximation  with 
unsteady,  incompressible  aerodynamic  forces,  incorporating  one  bending  mode 
and  one  torsion  mode,  are  shown  in  Fig  7.  The  tend  towards  coalescence  of  the 
bending  and  torsion  mode  can  be  seen  clearly  in  the  graph.  As  expected,  the 
torsional  mode  was  the  critical  flutter  mode,  which  crosses  the  imaginary  axis  at 
2 1 .44  m/s.  The  analytic  results  for  linear  flutter  theory  using  the  typical  section 
U-g  method  with  complex  sinusoidal  aerodynamic  force  formulation  is  shown  in 
Fig.  8.  Figure  9  shows  the  results  of  the  U-g  method  with  the  modifications  for 
varying  spanwise  deflection,  as  described  in  Section  2.4.  Again,  the  torsional 
mode  was  the  critical  flutter  mode  and  the  divergence  velocity  was  slightly  lower 
than  the  flutter  velocity. 

The  same  linear  results  can  also  be  obtained  by  using  the  harmonic  balance 
method  to  solve  the  complete  flutter  equations  formulated  in  Chapter  3.  This  was 
accomplished  by  inserting  the  root  angle  of  attack,  0r,  as  zero,  and  eliminating  all 
nonlinearity  in  the  formulation.  Comparison  of  results  from  the  various  analytic 
methods  shown  in  Table  1  are  seen  to  be  essentially  identical. 
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As  shown  in  Table  1,  the  linear  flutter  velocity,  Vp,  from  all  four  analysis 
methods  are  higher  than  the  linear  divergence  velocity,  Vp,  because  the  given 
properties  characterize  a  composite  wing  without  bending-torsion  structural 
coupling,  and  both  the  elastic  axis  and  the  C.G.  are  at  the  mid-chord.  The  trend 
towards  coalescence  of  the  bending  and  torsion  frequencies  can  also  be  seen 
clearly  in  Fig.  8  and  Fig.  9  -  dropping  of  the  torsional  frequency  from  about  24.2 
Hz  to  12.98  Hz  and  eventually  to  8  Hz,  while  the  bending  frequency  drops  from 
around  4.2  Hz  to  0  Hz  as  it  approaches  the  divergence  velocity. 

The  results  for  nonlinear,  large  amplitude  flutter  theory  using  the  harmonic 
balance  method  are  shown  in  Figs.  10  to  12.  The  figures  are  paired  in  groups  of 
two.  Fig.  10  shows  the  graphs  of  static  position  assuming  no  flutter,  and  the 
average  wing  angle  from  static  to  full  stall  flutter,  versus  velocity.  Both  the  static 
position  and  the  average  wing  angle  represent  the  sum  of  root  angle  of  attack  and 
the  twist  at  the  75%  span  location  due  to  aerodynamic  force.  Fig.  11  shows  the 
limit  cycle  oscillation  amplitude  versus  velocity  and  frequency  with  constant  root 
angle  of  attack  (0r  =  1°,  5°,  10°,  15°)  .  Fig.  12  shows  flutter  velocity  and 
frequency  versus  increasing  root  angle  of  attack. 

In  Fig.  10,  for  each  curve  of  constant  root  angle  of  attack,  both  the  steady, 
static  analysis  (solid  lines),  and  full,  unsteady  flutter  analysis  (dashed  lines)  are 
presented,  so  as  to  show  where  the  flutter  starts.  The  trend  shows  that  the  average 
wing  angle  exhibits  a  sharp  decrease  when  the  velocity  is  increased  past  the  flutter 
boundary  for  low  root  angle  of  attack,  0r=  1°  ,  while  for  larger  root  angle  of 
attack,  0r=5°,  10°,  and  15°,  the  drop  is  more  gradual.  The  reason  for  the  sudden 
drop  at  low  root  angle  of  attack  is  because,  when  the  root  angle  of  attack  is  below 
static  stall  angle  the  behavior  is  still  essentially  governed  by  the  linear 
aerodynamics,  so  the  oscillation  growth  is  still  essentially  exponential  up  until 
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deep  stall  is  reached,  therefore,  a  small  change  in  velocity,  would  result  in  large 
increase  in  total  angle.  The  characteristics  of  static  position  for  large  root  angle  of 
attack  are  essentially  the  same  as  for  the  low  root  angles  of  attack,  except  that  the 
change  from  light  stall  to  deep  stall  is  more  gradual,  since  the  wings  are  already  in 
deep  stall  from  the  initial  angle  of  attack  input. 

In  Fig.  11,  the  trend  shows  decreasing  flutter  velocity  with  increasing  root 
angle  of  attack.  Notice,  for  the  6r=1°  case,  the  velocity  began  to  decrease  at 
around  an  oscillation  amplitude  of  12°.  This  is  due  to  the  oscillation  penetrating 
the  negative  side  of  the  lift  deviation  curve.  The  same  trend  can  be  observed  from 
the  0r=5°  curve.  Also  notice,  for  root  angle  of  attack  of  10°  and  15°  at  low 
oscillation  amplitude,  the  velocity  is  constant  up  to  a  point.  This  is  because  the 
root  angle  of  attack  here  is  in  the  stall  range  initially,  and  while  the  oscillation 
amplitude  is  small,  the  oscillation  range  stays  only  in  the  stall  range.  When  the 
oscillation  become  large  again,  the  oscillating  range  exits  and  enters  the  stall 
range,  and  the  trend  of  increasing  velocity  with  increasing  amplitude  returns. 

Fig.  1 1  also  shows  the  limit  cycle  oscillation  frequency.  Here,  the  frequency 
increases  with  increasing  root  angle  of  attack.  As  the  oscillation  amplitude  is 
increased,  the  frequency  drops  slightly  until  the  amplitude  penetrates  the  negative 
deviation  curve,  and  then  the  frequency  increases. 

The  analytic  flutter  boundaries  are  presented  in  Fig.  12.  These  were  defined 
as  oscillations  for  amplitudes  less  than  1°.  Figure  12  shows  that  the  flutter 
velocity  starts  at  the  linear  flutter  velocity,  however,  due  to  the  close  proximity  of 
linear  flutter  speed  and  linear  divergence  velocity,  it  immediately  begins  to  exhibit 
nonlinear  behavior.  An  increase  in  the  root  angle  of  attack  causes  the  flutter 
velocity  to  drop  and  the  flutter  motion  to  become  more  purely  torsional.  This  can 
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also  be  seen  in  Fig.  12,  where  the  flutter  frequency  migrates  closer  to  the  natural 
torsion  frequency. 

The  results  and  trends  obtained  here  seem  to  correlate  generally  with  the 
experimental  results  obtained  by  Dunn  [Ref.  14].  Further  closer  correlation  would 
probably  require  improved  aerodynamic  fitting  of  the  appropriate  Reynold's 
number  aerodynamics,  and  inclusion  of  an  additional  structural  nonlinearity. 


Methods  k 

Vp  (m/sec) 

cOpfHz) 

Vofm/sec) 

Fade  Approx.* 

0.278 

21.44 

13.56 

20.83 

U-g  method* 

0.278 

21.43 

13.55 

20.82 

U-g  method** 

0.263 

21.69 

12.98 

20.82 

HBM** 

0.263 

21.70 

12.95 

20.86 

*  Typical  Section 

*♦  Varying  spanwise  deflection 


Table  1 .  Linear  flutter  characteristics 
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Figure  9.  U-g  flutter  analysis  (Varying  spanwise  deflection) 


Frequency  (Hz)  Oscillation  Amplitude  (deg) 
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CHAPTERS 


CONCLUSIONS 


The  present  investigation  has  developed  a  simple  aeroelastic  flutter  analysis 
using  the  ONERA  aerodynamic  model  to  predict  stall  flutter  behavior  of  a  3- 
dimensional  wing.  The  formulation  included  nonlinearities  in  aerodynamic 
effects,  and  uses  the  harmonic  balance  method  together  with  a  Newton-Raphson 
numerical  solution  technique. 

In  keeping  with  the  objective,  some  simplifications  were  made  for 
convenience,  but  the  model  still  retained  the  essential  nonlinear  characteristics  of 
the  aerodynamics.  For  example,  a  single  break  point  approximation  of 
aerodynamic  force  curves  was  used,  only  first  bending  and  first  torsion  mode 
shapes  were  considered,  the  angle  of  attack  at  75%  span  was  taken  as  the  average 
angle  of  attack  for  the  entire  wing,  and  only  first  harmonic  terms  in  the  harmonic 
balance  method  were  used.  These  modifications  can  be  extended  to  more 
complex  formulations  as  need  arises. 

As  shown  in  Chapter  4,  the  results  of  the  analytical  method  exhibit  almost 
all  trends  from  experimental  data  and  observation,  as  found  for  example  by  Dunn 
in  [Refs.  14  and  15].  The  trends  show  the  flutter  speed  decreases  as  the  root  angle 
of  attack  increases.  The  results  also  show  the  flutter  frequency  increases  toward 
the  torsional  natural  frequency  as  the  root  angle  of  attack  increases.  The  analysis 
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produced  nonlinear  limit  cycle  oscillations,  specifically,  it  showed  the 
"hardening"  trend  as  the  amplitude  of  limit  cycle  oscillation  increases. 

The  present  study  also  investigated  the  fitting  of  the  ONERA 
aerodynamic  model  to  measured  experimental  aerodynamic  stall  flutter 
characteristics.  The  results  described  in  Appendix  B  show  approximate 
representation  of  the  2-dimensional  nonlinear  aerodynamic  force  up  to 
certain  reduced  frequencies,  k,  and  oscillation  amplitude  Oy. 

The  current  investigation  concentrated  specifically  on  the  nonlinearity 
of  aerodynamic  forces.  A  structural  nonlinearity  can  also  be  readily 
incorporated  into  the  formulation  in  a  similar  manner,  by  this  harmonic 
balance  procedure,  as  discussed  briefly  in  Appendix  F.  This  would  further 
improve  correlation  between  analysis  and  experiment. 

It  should  be  mentioned  before  closing  that  for  simplicity,  the  analyses 
in  this  report  were  all  done  in  the  frequency  domain  using  the  harmonic 
balance  method.  As  such,  they  are  capable  of  capturing  the  steady  nonlinear 
limit  cycles  that  often  appear  at  larger  amplitudes  of  oscillation.  An 
alternative  approach  to  nonlinear  flutter  behavior  can  be  done  in  the  time 
domain,  using  the  same  ONERA  aerodynamic  equations  and  integrating  the 
equations  in  time.  See  Tang  and  Dowell  [Ref.  31,  32].  These  time  domain 
solutions,  when  integrated  for  long  enough  times,  often  do  settle  down  to  the 
steady  nonlinear  limit  cycle  solutions  investigated  in  this  report.  However, 
some  nonlinear  regimes  may  also  lead  to  unbounded  solutions,  while  other 
regimes  may  lead  to  boimded,  but  nonperiodic,  chaotic-type  solutions.  These 
chaotic  solutions  can  occur  for  example,  when  oscillations  get  trapped  at  the 
end  points  of  a  free-play  structural  nonlinearity,  and  are  very  sensitive  to  the 
initial  starting  conditions.  See,  for  example^  Tang  and  Dowell  [Refs.  31,  32], 
Hauenstein,  Zara,  Eversman,  and  Qumei  [Ref.  33]  and  Yang  and  Zhao  [Ref.  34] 
for  more  details  on  these  nonlinear  solutions. 
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APPENDIX  A 

COEFFICIENTS  OF  AERODYNAMIC  FORCES 


This  section  lists  the  coefficients  of  the  aerodynamic  force  models  used  in 
present  study.  Some  of  the  coefficients  were  derived  from  aerodynamic  theory, 
others  were  obtained  from  previous  works  [Ref  14]. 

The  linear  portion  of  ONERA  ,  2-dimensional  aerodynamic  formulation,  are 
described  by  equations  (3.2)  and  (3.3).  The  linear  coefficients  (s^j,  Sz2>  Sz3,  ?ii,  A.2, 
and  a^z)  were  derived  from  2-dimensional  incompressible  flow  theory.  In  the 
complete  aerodynamic,  incompressible  flow  theory,  lift  can  be  formulated  as, 


I  =  +  iye-dae] 


+ 


2npCrdC(Ji:) 


- a 

-h  +  UQ  +  b 

e 

I2  j 

(A.1) 


where  a  is  a  nondimensional  parameter  defining  the  distance  from  the  mid-chord 
to  the  elastic  axis  divided  by  half-chord  length  (aft  direction  is  defined  as 
positive)  For  the  present  study,  a  =  -1/2  and  the  h  and  0  above  are  consequently 
defined  at  the  quarter-chord.  Again,  define  the  non-dimensional  time  variable  x, 
where  the  non-dimensional  time  derivative  is  denoted  by  ”  *  Equation  (A.l) 
then  becomes. 
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with  some  algebraic  manipulation  equation  (A.2)  becomes, 


i  =  ^pt/=c(Q,-Q,)  +  ipC/=cC„ 
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where  Cl  i-Cly  and  Cz^  are, 
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Substituting  value  of  a  and  equation  (2.5)  into  equation  (A.4)  gives. 


•  ** 

(A.6) 

where  the  angle  of  attack  a  =  Q  -  h  /  b.  Place  the  single  lag  pole  approximation 
of  the  Theodorsen  fiinction,  defined  in  Eqs.  (3.25)  and  (3.26),  into  equation  (A.5). 
Similarly,  substitute  equation  (2.5)  into  equation  (A.5).  This  results  in  the 
following  equation  fi-om  which  the  coefficients  can  then  be  found, 
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The  coefficients  of  linear  moment  formulations  were  found  in  a  similar  way.  The 
linear  coefficients  are  summarized  in  Table  2. 

The  values  of  aerodynamic  integrals  appeared  in  equations  (3.134),  (3.135) 
and  (3.144)  to  (3.148)  are  shown  in  Table  3.  The  integrands  are  the  mode  shapes, 
or  combinations  of  mode  shapes  for  a  cantilevered  beam.[Ref  26  and  27  ].  Note 
that  for  simplicity,  only  first  bending  and  torsion  modes  were  used. 


h  =  Jo  <t>/,  <l>a  ^ 

^  =  Id  <t>a 

h  =  Jd  ^ 

h  =  Jd  <l>a  ^ 


(A.8) 

(A.9) 

(A.10) 

(A.11) 

(A.12) 


Additionally,  as  mentioned  in  Chapter  3,  for  the  nonlinear  portion  of  the 
aerodynamic  force,  stalling  was  only  evaluated  at  the  75%  span  point.  It  was 
decided  that  the  behavior  at  75%  span  is  typical  of  the  entire  wing.  This  idea 
relates  back  to  the  2-dimensional  typical  section  analysis  described  in  Section  2.1. 

The  material  properties  and  characteristics  were  obtained  from  previous 
study.  These  properties  describes  the  characteristics  of  the  [03/90]s  graphite/epoxy 
test  specimen  (see  table  4). 


Table  2.  Linear  aerodynamic  coefficients  (a  =  -7/2) 


1.0000  I  0.6779  I  0.5000  |  0.7830  |  0.6366 

Tables.  Aerodynamic  integrals 


M  =  0.283  kg/m 

p  =  1 .23  kg/m2 

S„  =  0.0kg 

c  =  0.140  m 

=  0.343  e  -3  kg-m 

I  =  0.559  m 

Kh  =  206.4  N/m2 

b  =  0.070  m 

Ke,=  8.209  N/rad 

©h  =  27.02  rad/sec 

S  =  0.783  m2 

©a  =  154.6  rad/sec 

Table  4.  Specimen  properties  from  previous  study 


APPENDIX  B 

FITTING  OF  ONERA  AERODYNAMIC  MODEL 


Fitting  of  the  2-dimensional  ONERA  aerodynamic  coefficients  to  the 
experimental  lift  and  moment  data  of  the  oscillating,  stalled  NACA  0012  airfoil 
given  in  Ref.  24  is  explained  in  this  section.  This  is  done,  first,  to  determine  the 
coefficients  which  correctly  model  the  aerodynamic  forces.  Second,  the  fitting  of 
the  aerodynamics  allows  for  closer  examination  of  the  effects  of  each  of  the 
harmonics  on  the  aerodynamic  hysteresis  cycle. 

The  aerodynamic  force  coefficients  can  be  written  in  harmonic  form  by 
assuming  harmonic  motion,  as  described  in  Section  3.3.  The  nonlinear 
coefficients  in  equation  (3.4),  r„  rj  and  Tj,  were  determined  by  Petot  and  Loiseau 
[Ref  16]  for  an  OA  209  airfoil.  These  coefficients  were  later  modified  by  Dunn 
and  Dugundji  [Ref  15].  All  of  the  previous  formulation  of  these  coefficients 
were  extremely  complex.  In  this  study,  adhering  to  the  objective  of  simplifying 
nonlinear  analytical  methods,  these  coefficients  were  simplified.  They  are  of  the 
form. 


fi  = 


^2  = 


^3  ~  ^30  ^32  ( AG20  ) 


(B.l) 

(B.2) 

(B.3) 
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The  major  simplification  comes  in  terms  of  using  only  the  constant  term  of  the 
nonlinear  deviation  factor,  AC^,  where  previous  works  had  used  the  full  terms. 
The  coefficients,  rjo,  rj2,  r2o,  r22,  r3o,  and  r32,  are  then  determined  by  fitting  to  the 
experimental  data  by  McAlister,  Pucci,  McCroskey  and  Carr  [Ref  24].  The 
coefficients  along  with  the  stall  angle,  aj,  and  the  slope  of  nonlinear  deviation,  b„ 
are  summarized  in  Table  5.  Table  6  shows  the  reduced  frequency,  k,  the  vibratory 
angle,  a^,  and  the  angle  of  attack,  a^,  for  each  of  the  experimental  cases  to  be 
fitted.  Note,  when  transferring  aerodynamic  forces  (Section  3.4b)  from  quarter 
chord  to  mid  chord,  the  effect  of  moment  is  small  compared  to  that  of  the  lift. 
Therefore,  only  the  lift  is  fitted  to  the  experimental  results.  In  the  present 
analysis,  the  coefficients,  rj,  used  for  the  moment  is  the  same  as  those  used  for  the 
lift. 


Fourier  analysis  was  performed  to  extract  the  first  and  second  harmonics  of 
the  experimental  data  in  preparation  for  fitting  comparison  (see  Table  7).  The 
resulting  aerodynamic  hysteresis  loop  from  the  fitting  process  are  plotted  against 
the  experimental  data  in  Fig.  13.  In  addition,  the  static  approximation  of  lift  and 
moment  curve  used  in  the  fitting  process  is  shown  in  Fig.  14. 

Each  individual  harmonic  component  affects  the  appearance  of  the  hysteresis 
cycle  distinctly.  The  gross  characteristics  of  each  harmonic  are  shown  in  Fig.  15. 
Fig.  15a  shows  the  effect  of  the  first  sine  harmonic  to  the  hysteresis.  It  generally 
follows  along  the  linear  force  curve.  The  first  cosine  harmonic  influences  the 
hysteresis  cycle  on  the  amount  of  deviation  from  the  linear  curve  due  to  static 
stalling,  as  shown  by  Fig.  15b.  If  the  second  harmonic  sine  term  is  large  enough, 
a  crossing  of  the  hysteresis  loop  appears  as  shown  by  Fig.  15c.  Finally,  the 
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second  cosine  term  determines  the  degree  of  curvature  of  the  resultant  hysteresis 
loop  (see  Fig.  15d).  Collectively,  each  of  these  harmonic  terms  influence  the 
cycle.  However,  these  are  just  broad  generalities.  In  reality,  other  physical 
aspects  also  influence  the  nature  of  fluid  flow.  Nevertheless,  this  exercise  gives 
insight  into  the  basic  building  blocks  of  the  appearance  of  the  hysteresis  loops. 


r,o  =  0.700  ri2  =  0.150  a,  =10° 

r2o  =  0.246  r22  =  0.005  b,  -  lift  =  8.30 

1*30  ~  -0.024  r32  =  0. 1 1 6  b,  -  moment  =  0.40 

T able  5 .  Aerodynamic  fitting  parameters 


Case 

“o 

tty 

k 

1 

9.8° 

o 

o^ 

0.104 

2 

9.8° 

9.9° 

0.151 

3 

bo 

o 

9.9° 

0.253 

4 

14.9° 

9.9° 

0.104 

5 

14.9° 

VO 

o 

0.153 

6 

14.7° 

14.0° 

0.104 

Table  6.  Aerodynamic  fitting  specifications 
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Experimental  Values  — 


Case 

s, 

Cl 

C2 

1 

0.818 

0.747 

0.420 

0.159 

-0.123 

2 

0.883 

0.748 

0.347 

0.072 

0.022 

3 

0.983 

0.801 

0.389 

0.113 

-0.073 

4 

0.897 

0.497 

0.574 

0.013 

-0.027 

5 

0.868 

0.662 

0.519 

-0.014 

0.063 

6 

0.826 

0.661 

0.506 

0.056 

0.016 

Analytical  Results  — 


Case 

Co 

s, 

c. 

s, 

Cj 

1 

0.632 

0.688 

0.316 

0.121 

0.012 

2 

0.632 

0.751 

0.312 

0.083 

-0.012 

3 

0.632 

0.766 

0.338 

0.039 

-0.021 

4 

0.764 

0.503 

0.504 

0.071 

0.024 

5 

0.764 

0.620 

0.468 

0.051 

0.009 

6 

0.588 

0.771 

0.581 

0.122 

0.055 

Table  7.  Harmonic  components  of 
aerodynamic  fitting  results 
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Experimental  Values  - 


Case 

n 

o 

s, 

c, 

§2 

C2 

1 

-0.090 

-0.097 

-0.017 

0.016 

0.080 

2 

-0.125 

-0.103 

-0.002 

0.038 

0.063 

3 

-0.128 

-0.074 

0.015 

0.080 

0.013 

4 

-0.058 

-0.130 

-0.036 

-0.025 

0.063 

5 

-0.037 

-0.142 

-0.017 

0.003 

0.051 

6 

-0.086 

-0.131 

-0.042 

-0.02 

0.065 

Table  7  continued,  Harmonic  components 
of  aerodynamic  fitting  results 


Coefficient  of  lift  Coefficient  of  lift 
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a.  Case  1 


b.  Case  2 


Figure  13.  Aerodynamic  fitting 


Coefncient  of  lift  Coefficient  of  lift 
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Angle  of  attack  (deg) 


c.  Case  3 


d.  Case  4 


Figure  13.  continued,  Aerodynamic  fitting 


Coefficient  of  lift  Coefficient  of  lift 
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Analysis 


Experiment 


“„=  14.9° 

“  =9.9° 
k  =  .153 


0  5  10  15  20  25  30 


Angle  of  attack  (deg) 


e.  Case  5 


f.  Case  6 


Figure  13.  continued,  Aerodynamic  fitting 


Coefficient  of  moment  Coefficient  of  lift 
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a.  Static  lift  curve  approximation 


b.  Static  moment  curve  approximation 


Figure  14.  Static  aerodynamic  force  curves 
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b.  First  harmonic  cosine  term 


c.  Second  harmonic  sine  term 


F igure  1 5 .  Effects  of  harmonics 


APPENDIX  C 

EXPERIMENTAL  DATA  FOR  FITTING 


As  a  reference  for  the  fitting  study  described  in  the  previous  section,  the 
original  experimental  data  from  McAlister,  Pucci,  McCroskey,  and  Carr  [Ref.  24] 
is  included  here. 
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0012  AiRrOlt 


Figure  12.-  Contiitw'ed 


NACA  0012  AiRFOlt 
PftAME  :  6102 


Figure  12.-  Continued 


APPENDIX  D 

AERODYNAMIC  FORCE  CURVES 

Equations  (3.41)  to  (3.43)  gives  the  nonlinear  aerodynamic  force  deviation 
curve,  where, 

Ti  [{  a,.  2  J 

AC,.=^[-U,-^l-jsin(2(l).)l  (D.2) 

7t  L  V  ^  ^ 

ACzc2  =  cos(t),  - Sin3<|),  (D.3) 

71  L  2  6  J 

The  non-dimensional  time  parameter,  (j)],  associated  with  the  above 
formulation  are  still  applicable,  as  stated  in  Eqs.  (3.44)  and  (3.45), 


sin'^  6 

if 

-1<6<1 

71 

2 

if 

8>  1 

(D.4) 

71 

if 

5  <  -1 

2 


where. 
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(D.5) 


This  single  break  point,  simplified  model  of  the  lift  curve  is  shown  in  Figure  3.  In 
the  present  study,  with  the  combination  of  root  angle  of  attack  and  oscillation 


amplitude,  the  range  of  oscillation  amplitude  may  reach  the  negative  portion  of 
the  force  curve.  As  a  result,  symmetric  aerodynamic  force  curves  were  used.  The 


negative  portion  of  the  force  curves  were  accounted  for  by  including  a  second 
stall  angle  d,,  where  d,  =  -a,.  The  negative  portion  of  the  aerodynamic  force 


was  then  combined  with  the  single  sided  formulation  to  yield  the  expanded 
version  of  Eqs.  (D.l),  (D.2)  and  (D.3), 
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where  the  addendum  to  Eqs.  (3.34)  and  (3.35)  are, 

if  -  1  <  8  <  1 

if  8  >  1  (D.9) 

if  8<-l 

ctr 

Fig.  16  shows  the  resulting  symmetric  aerodynamic  force  curve. 


(D.IO) 
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Figure  16.  Symmetric  aerodynamic  force  curve 


APPENDIX  E 

NEWTON-RAPHSON  METHOD 


The  Newton-Raphson  method  is  a  powerful  iterative  algorithm  for  solving 
nonlinear  equations.  This  method  gives  one  a  very  efficient  means  of  converging 
to  a  solution,  provided  that  a  sufficiently  good  initial  guess  has  been  established. 

A  typical  problem,  deals  with  N  functional  equations  to  be  zeroed,  involving 
variables  Xj,  where  i=l,2,;..,  N, 

=  0  (E.l) 

where  x  denotes  the  state  vector  and  F  denotes  the  entire  vector  of  residual 
functions.  The  functions  F  can  be  expanded  in  Taylor  series 


F;(x  +  6x)  =  F;(x)+ +o(6x')  (E.2) 

j=idXj 

The  matrix  of  partial  derivatives  in  above  equation  is  the  Jacobian  matrix  J.  By 
neglecting  terms  of  order  6jc^  and  higher  and  setting  F(x  +  5x)  =  0,  the  resulting 
equations  become  a  set  of  linear  equations  which  define  the  correction  factor  5x 
that  move  each  function  closer  to  zero  simultaneously, 

8x,.,  =  -y-'F(x,„)  (E.3) 


The  corrections  are  then  added  to  the  solution  vector, 

^(n+l)  ~  ^(n)  ^^(n) 
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(E.4) 
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and  the  process  is  iterated  to  convergence. 

The  Newton-Raphson  iteration  method  is  applied  to  the  present  analysis  by 
rearranging  equation  (3.89)  as 

[M]{q} +  [K]{q} -{Q}  =  0  (E_5) 

Equation  (E.5)  is  the  matrix  form  of  a  system  of  six  nonlinear  equations  with  six 
unknowns.  The  unknowns,  or  the  state  vector  x  is  comprised  of  the  harmonic 
components  of  the  modal  amplitudes,  qj^,  qj^,  and  qj^..  Some  adjustment  of  the 
state  vector  is  made  to  ensure  convergence  to  the  non-trivial  solution.  The 
adjustment  involves  setting  the  sine  component  of  one  mode,  q^g  to  a  small 
constant,  which  defines  the  amplitude  level  of  the  oscillation,  while  its  cosine 
component,  q2c,  is  set  to  zero,  since  the  flutter  limit  cycle  oscillations  can  start  at 
any  arbitrary  phase.  The  mode  chosen  was  the  torsional  mode,  since  the  torsional 
mode  dominates  the  oscillation  behavior.  These  two  components  were  replaced 
by  ka  and  k,  the  reduced  torsional  frequency  and  reduced  flutter  frequency 
respectively.  The  resulting  state  vector  x  is  then. 


^10 

^20 

Qxs 

Qxc 

K 

k 


(E.6) 
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The  Jacobian  matrix  was  solved  numerically.  This  was  done  by  increasing 
the  state  vector  x  a  small  amount,  finding  the  incremental  change  in  the  residual 
vector  F,  and  then  evaluating  the  partial  derivative  dVJ  dx^-,  as  AFj/  Axj. 

The  Newton-Raphson  procedure  for  solving  flutter  problems  by  introducing 
k  and  into  the  state  vector  x  described  here,  is  similar  to  that  used  by  Kuo, 
Morino,  and  Dugundji  [Ref.  29]  for  nonlinear  panel  flutter. 


APPENDIX  F 

INCLUDING  STRUCTURAL  NONLINEARITIES 


The  inclusion  of  structural  nonlinearities  is  relatively  straightforward 
with  the  Harmonic  Balance  Method.  An  example  of  a  simple  cubic  torsional 
nonlinearity  was  included  in  the  analysis  of  Dunn  &  Dugundji  [Ref.  15].  The 
present  appendix  shows  a  simple  cubic  torsional  structural  nonlinearity 
included  in  the  previous  analyses  of  this  report.  Other  types  of  structural 
nonlinearities  can  also  be  simply  incorporated. 

Turning  to  the  nonlinear  flutter  analysis  of  Section  3.4,  one  modifies 
the  basic  equations  of  motion  Eqs.  (3.89)  to  read. 


q,  - 

(F.l) 

M22  q2  +  K22  (i  +  Cn  e|j)  q2  =  Q2 

(F.2) 

where  the  structural  twist  angle  at  the  tip,  GgT/  is  related  to  the  twist 
coordinate  qi,  through  the  relation, 

Ssr  =  il'amq2  =  K 

See  Eqs.  (3.78),  (3.80),  and  (3.85).  The  nonlinear  stiffness  constant,  Cn,  is  of  the 
order  of  20  for  die  experimental  parameters  given  in  Table  4.  After 
nondimensionalization,  Eqs  (F.l)  and  (F.2)  lead  to  the  nonlinear  structural 
counterparts  of  Eqs.  (3.103)  and  (3.104),  namely. 
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^7Cli(qi  +  flPki 


:kJ=f 

Jo 


^Li/2  dx 


(F.4) 


^ith[^2+  ^092)  =  j*^  <l>a<^Mi/2dx  - 


i^IskaCNl^ 


(F.5) 


Introducing  Harmonic  Balance, 


Qi  =  qio  +  qissinkx  +  ^ccoskx 


(F.6) 


leads  to  the  additional  nonlinear  terms  on  the  right  hand  sides  of  Eqs.  (3.145), 
(3.148)  and  (3.149)  respectively,  namely. 


-  HAlil3ktCN(qi  +  q2c])q2o 

(F.7) 

-  n-|^4i3iq.CN(3  4,  +  |[q2s  +  4J)q2s 

(F.8) 

-^l^■iI3l4cN(3qi+  |[te+q2e])q2e 

(F.9) 

The  new  modified  nonlinear  Eqs.  (3.144)  to  (3.149)  can  then  be  solved  using 
the  Newton-Raphson  iteration  method,  as  described  previously. 

For  a  more  general  structural  nonlinearity,  Fn  =  Fn  (a,  a),  one  can 
characterize  the  nonlinearity  as  in  the  previous  aerodynamic  case,  by 
assuming  the  motion  a  to  be  described  as 
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a  =  Oo  +  otv  sin  <))  (F.IO) 

See  previous  Eq.  (3.33).  Then  by  placing  this  motion  back  into  Fn,  one  can 
represent  the  corresponding  nonlinear  force  as, 

Ffj  =  ^No  "*■  ®NS1  ^  ®NC1  COS<j)  (F.ll) 

where  the  Fourier  coeffecients  are  given  as. 


B 


('In 

Fn  (a, a)  d(t) 

0 

rln 


®NS1 


^NCl 


JFivj(a,a)  sin(})  dcj) 
0 

J'2n 

Ffj  (a,a)  cos  <j)  d<t) 

0 


(F.12) 


(F.13) 


(F.14) 


This  is  analogous  to  Eqs.  (3.37)  to  (3.40).  Finally,  upon  noting  that 
<|)  =  kr  +  ^ ,  one  can  revert  back  to  the  basic  forms, 

a  =  +  as  sin  kx  +  Oc  cos  kx  (F.15) 

,  =  Fjmo  +  Fj^si  sin  kx  +  Fnci  cos  kx  (F.16) 

as  in  the  previous  Eqs.  (3.67)  to  (3.71).  These  nonlinear  forces  Eq.  (F.16)  are 
then  easily  incorporated  into  the  Harmonic  Balance  analysis. 


APPENDIX  G 
ANALYSIS  CODES 


In  the  process  to  complete  this  study,  many  programs  were  written.  A  small 
portion  of  them  were  written  as  a  learning  tool  for  the  author,  others  are  building 
blocks  for  the  final  full  flutter  analysis.  Furthermore,  in  order  to  save  precious 
time,  some  codes  were  compiled  several  times  with  only  minute  differences. 
Thus,  only  those  codes  which  are  essential  to  aiding  the  completeness  of  this  thesis 
are  listed  here. 

All  codes  were  in  FORTRAN  and  compiled  with  Microsoft  FORTRAN™ 
compiler  Ver.  5.10  on  an  IBM  compatible  386,  25  MHz  personal  computer.  Some 
subroutines  were  taken  fi-om  Numerical  Recipe  book  [Ref  28].  They  are  also 
available  on  Athena  computing  system. 
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cccccccccc 


c  AEROLl.FOR  Aerodynamic  forces.  Harmonic  Balance 
c  method  used, 

c  -  Linear  Lift,  CLl 

c  Warren  Chen,  2-19-93 
cccccccccc 

parameter  (nn=2) 
double  precision  cll 
real  phi, k,phid, theta 
k;=:0.2 

open { 5 , f ile= ' aerol .in', status= ' old ' ) 

open(9,file='aeroll.out' , status= ' unknown ' ) 

write ( 9 , * ) '  AERODYNAMIC  FORCES ' 

write(9,*)'  LI FT -LINEAR  PORTION' 

write (9, 1000) 'K=  ',k 

write  (9, *)  '  ' 

write (9, *) '  ' 

write (9, *) 'THETA{deg)  Cll' 

write (9,*) ' - • 

10  continue 

read (5, *)phid 
if  (phid.ge.999.9)  then 
goto  500 
endif 

phi= (phid*3 .1416) /ISO . 0 
call  funcv (k, phi, cll, theta) 
thetad=theta*180 . 0/3 . 14159 
write (9,2000) thetad, '  ',cll 

goto  1,0 
500  continue 

write (9,*) 'End  of  Data' 
stop 

c  formats 

1000  formate  ',a,f7.4) 

2000  formate  ' ,  f5 .  l,  a,  f6 . 4) 

3000  formate  ',a,f8.6) 
end 


c  Subroutine  to  calculate  lift 
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subroutine  f uncv(k, phi, p, theta) 

real  sll , sl2 , sl3 , aol , laml , lam2 , tho, ths , the , 
$  phi, k, he, hs,b, theta 
real*8  p, Is, le, f , g, elgo, elgs, clge 


e  Variable  initialization 


b=.07 
sll=2.142 
Sl2=1.571 
Sl3=0.0 
aol=5.9 
laml=0 . 15 
lam2=0 . 55 
tho=0 . 0 
ths=0 .175 
the=0 . 0 
he=0 . 0 
hs  =  0 . 0 

ls=aol* (ths-k*he/b-k*the) 

le=aol* (the+k*hs/b+k*ths) 

f = (k**2*lam2+laml**2) / (laml**2+k**2) 

g=k*laml* (lam2-l) / (laml**2+k**2) 

elgo=aol*tho 

elgs=f *ls-g*le 

elge=g*ls+f *le 

theta=tho+ths*sin (phi) +the*eos (phi) 
p= (sll* (-k*thc- (hs/b) *k**2) -sl2*ths*k**2 
$  -sl3*k*the+elgs) *sin (phi) , 

$  + (sll* (ths*k- (he/b) *k**2 ) 

$  -sl2*thc*k**2+sl3*k*ths+elge) *eos (phi) 

$  +elgo 

return 

end 


cccccccccc 

c  AER0L4.F0R  Aerodynamic  forces.  Harmonic  Balance 
c  method  used, 

c  ~  Total  Lift,  C11+C12 

c  New  Basic  Case.  Taken  from  McAlister's  NASA 
c  experimental  work. 

c  New  parameters (stall  angles,  lift  curve  slopes, 
c  more  accurate  rl,r2,r3's 

c  by  taking  into  account  variability  of  Dcz's 

c  Trying  to  get  a  better  fitting  by  vary 

c  rl,r2,r3. 

c 

c  Warren  Chen,  3-7-93 

cccccccccc 

c 

real*8  phi , k, phid, theta, cz2c2 , cz2s2 ,  f ,  g 
$  , cz2cl, cz2sl, cz20, clgo, els, clc, cll, cl2 
$  , alphaOd, alphavd, alphaO , avib, rl , r2 , r3 
write (* , *) ' cccccccccccccccccccccccccccccccccccc 
write (*,*)' cThis  program  calculates  the  first  c 
write (*,*)' cand  second  harmonic  components  of  c 
write (*,*) 'cthe  2-d  liftloop.  c 

write (*,*) 'c  The  program  inputs  are:  c 

write (*,*) 'c  reduced  frequency,  k,  c 

write(*,*)'c  average  wing  angle, in  degrees  c 

write(*,*)'c  vibration  angle, in  degrees  c 

write{*,*)'c  The  output  is  saved  in  file  c 

write(*,*)'c  AER0L4.0UT  c 

write (*,*) ' cccccccccccccccccccccccccccccccccccc 
write (*,*)'' 
write (*,*)' Enter  k: ' 
read{*,  *)k 

write (*,*) 'Enter  alphaO  in  deg:' 

read (*,*) alphaOd 

write (*,*)' Enter  alphav  in  deg:' 

■  read ( * , * ) alphavd 

alpha0=alpha0d*3 . 1415926/180 . 0 

avib=alphavd*3 . 1415926/180 . 0 

open ( 9 , f ile= ' aerol4 . out ' , status= ' unknown ' ) 

write (9,*)'  AERODYNAMIC  FORCES - TOTAL  LIFT 

$  AEROL4 . out ' 

write (9, 1000) ' INPUT:  K=  ',k,'alpha0=  ' 
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$  ,alphaOd, '  alphav=  ',alphavd 

write (9, *) '  ' 

call  funcvl (alphaO , avib, k, phi, theta, clgo, els 
$  ,elc,f ,g,cll) 

call  f unevn ( alphaO , avib , rl , r2 , r3 , k , phi 
$  , cz2c2 , cz2s2 , cz2cl, cz2sl , cz20 , cl2) 
write(9,5000) 'F=  ',f,'  G=  ',g 

write(9,6000) 'rl=  ',rl,'  r2=  ',r2,'  r3= 

$  '  ,r3 

write  (9 , *)  '  ' 

write (9, *) 'Linear  portion,  Cll' 

write  (9, 4000)  ' cllO=  ',clgo,'  clls=  ',cls,' 

$  cllc=  ' , clc 
write(9,*)'  ' 

write  (9, *)  'NonLinear  portion,  C12 ' 

write  (9 , 4000)  ' cl20=  ',cz20,'  cl2sl=  ',cz2sl,' 

$  cl2cl=  ' , cz2cl 

write (9, 5000) ' cl2s2=  ',cz2s2,'  cl2c2=  ',cz2c2 

write (9, *) '  ' 

write (9 ,*)' Complete  Cl' 

write(9,4000) 'c0=  ' , clgo+cz20 , '  sl=  ',cls+cz2sl 
$  ,'  cl=  ',clc+cz2cl 
write (9, 5000) ' s2=  ',cz2s2,'  c2=  ',cz2c2 

write (9 , *) '  ' 

write (9, *) 'PHI  THETA  CLl  CL2 

$  CL' 

write (9,*) ' (deg)  (deg) ' 

write  (9,*)' - ' 

c 

c 

phid=0 . 0 
10  continue 

if  (phid.ge . 361 . 0 )  then 

goto  500 

endif 

phi=(phid*3. 1415826) /180.0 

call  funcvl (alphaO , avib, k, phi, theta, clgo, els 
$  , clc, f ,g,cll) 

call  f unevn ( alphaO , avib , rl , r2 , r3 , k , phi , cz2c2 , 

$  cz2s2,cz2cl,cz2sl,cz20,cl2) 
thetad=theta*180 . 0/3 . 1415926 
write (9 , 2000) phid, thetad, cll , cl2 , cll+cl2 
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phid  =phid  +  15.0 
goto  10 
500  continue 

write (9, *) '  ' 
write (9, *) 'End  of  Data' 
stop 
c 

c  formats 


1000  format ( 
2000  format { 
3000  format ( 
4000  format { 
5000  format ( 
6000  format ( 
end 


' ,a,f6.3) 

'  ,f6.1,f8.3,f9.3,f9.3,f8.3) 
'  ,a,f8.6) 

'  ,a,f8.4,a,f8.4,a,f6.4) 

'  ,a,f6.4,a,f6.4) 

'  ,a,f5.3,a,f5.3,a,f6.3) 


c 

cccc 

c  Subroutine  to  calculate  the  linear  portion 

cccc 

c 

subroutine  funcvl (alphaO , avib, k, phi, theta, clgo 
$  ,pls,plc, f ,g,p) 

real *8  sll,sl2,sl3, aol , laml , lam2 , alphaO 
$  , avib , the , p , phi , k , he , hs , b , theta , pis , pic 
$  , Is, Ic, f ,g, clgo, clgs, clgc 
c 

b=..07 
sll=3 .142 
Sl2=1.571 
sl3=0. 0 
aol=6 .28 
laml=0 . 15 
lam2=0 . 55 
thc=0 . 0 
hc=0 . 0 
hs=0 . 0 
c 

ls=aol* (avib-k*hc/b-k*thc) 
lc=aol* (thc+k*hs/b+k*avib) 
f= (k**2*lam2+laml**2) / (laml**2+k**2) 
g=k*laml* (lam2-l) / (laml**2+k**2) 


127 


clgo=aol *alphaO 
clgs=f *ls-g*lc 
clgc=g*ls+f*lc 

theta=alphaO+avib*sin (phi) +thc*cos (phi) 
pls=  (sll*  (-}c*thc-  (hs/b)  *lc**2)  -sl2*avib*}c**2 - 
$  sl3*k*thc+clgs) 
plc=  (sll*  (avib*lc-  (hc/b)  *]c**2)  - 
$  sl2*thc*k**2+sl3*k*avib+clgc) 
p=pls*sin (phi) +plc*cos (phi) +clgo 
return 
end 
c 

cccc 

c  Subroutine  to  calculate  the  Nonlinear  portion 

cccc 

c 

subroutine  f uncvn ( alphaO , avib , r 1 , r2 , r3 , k , phi 
$  ,cz2c2, cz2s2,cz2cl,cz2sl,cz20,p2) 
double  precision  dczO , dczsl , dczc2 , rO , rsl , rcl , rs2 
$  ,kl,k2,k3,k4,cz2c2,cz2s2,cz2cl,cz2sl,cz20 
$  , alphal , bl , alphaO , avib , rl , r2 , r3 , phil , k , p2 , phi 
$  , pent2 , philt , alphalt 
c 

alphal=0 . 1745 
alphalt=-0 . 1745 
bl=8.3 
c 

pent= (alphal -alphaO) /avib 
if  (pent  .gt.  1.0)  then 
phil=1.571 

else  if  (pent  .It.  -1.0)  then 

phil=-l . 571 

else 

phil=asin (pent ) 
endif 
c 

pent2= (alphalt-alphaO) /avib 
if  (pent 2  .gt.  1.0)  then 
philt=1.571 

else  if  (pent2  .It.  -1.0)  then 

philt=-1.571 

else 


philt=asin (pent2 ) 
endif 

dczO= (bl*avib/3 . 14159) * (-pent* (1.571-phil) 
$  +COS (phil) ) - (bl*avib/3. 14159) * {pent2* 

$  (philt+1 . 571) +COS (philt) ) 
dc2sl= (bl*avib/3 .14159) * ( (1.571 -phil) - 
$  0 . 5*sin (2 . 0*phil) ) - (bl*avib/3 . 14159) * 

$  ( - (philt+1 .571) -0.5*sin(2 . 0*philt ) ) 
dczc2= (bl*avib/3 . 14159) * (-0 . 5*cos (phil) - 
$  0 . 16G667*cos (3 . 0*phii) ) 

$  - (bl*avib/3 .14159) * (-0.5*cos (philt) - 
$  0 . 166667*cos (3 . 0*philt ) ) 
rl=0 . 7+0 . 15*dcz0**2 
r2=  (0 . 2465  +  0 . 005*dcz0**2) **2 
r3  =  (-0 .4  +  1 . 91*dcz0**2) *(0.2465**2) 
kl=r2-]c*k 
k2=rl*k 
k3=r2-4*k*k 
k4=2*rl*k 
r0=-r2*dcz0 
rsl=-r2*dczsl 
rcl=-r3*k*dczsl 
rs2=2*r3*k*dczc2 
rc2=-r2*dczc2 
cz20=r0/r2 

cz2sl= (kl*rsl+k2*rcl) / (kl**2+k2**2) 
cz2cl= (kl*rcl-k2*rsl) / (kl**2+k2**2) 

C22s2= (k3*rs2+k4*rc2) / (k3**2+k4**2) 
cz2c2=<k3*rc2-k4*rs2) / (k3**2+k4**2) 
p2=cz20+cz2sl*sin (phi) +cz2cl*cos (phi) 

$  +C22s2*sin (2 . 0*phi) +cz2c2*cos (2 . 0*phi) 
return 
end 


ccccc 

c  FOURI.FOR'  Program  to  calculate  Fourier  Coefficients 
c  Output : 1st  and  2nd  harmonics  of  aero 

c  hysteresis  curve  values 

c  Warren  Chen  3-16-93 
ccccc 


parameter (n=16,m=8) 
integer  i 

real*8  f(ltn) |C0/Cl/Sl/C2^s2 
open (5, f ile= ' fouri . in' , status= ' old' ) 
open(9, file= ' fouri.out ' ,status=' unknown' ) 
do  i=l,n 
read (5 , *) f  (i) 
enddo 
c0=0.0 
cl=0 . 0 
c2=0.0 
sl  =  0 . 0 
s2=0 . 0 
do  i=l,n 
cO=cO+f (i) 

sl=sl+f (i) *dsin(6.2832* (i-1) /n) 
cl=cl+f (i) *dcos (6.2832* (i-l) /n) 
s2=s2+f (i)*dsin(6.2832*(i-l)/m) 
c2=c2+f (i) *dcos (6.2832* (i-l) /m) 
enddo 
cO=cO/n 
sl=sl/m 
cl=cl/m 
s2=s2/m 
c2=c2/m 

write (9, *) 'Fourier  Coefficients' 
write ( 9 , * ) ' Aerodynamic  Moments ' 
write(9,*)'  cO  si  cl  s2 

write (9, 1000) cO, sl,cl, s2,c2 
1000  format ( '  ' , f 7 . 3 , f 7 . 3 , f 7 . 3 , f 7 . 3 , f 7 . 3) 

end 
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cccccccccc 

c  HARMVdl.FOR  Complete  3-D  flutter  program, 
c  xfl)=ql0,  x(2)=q20,  x(3)=qls,  x(4)=qlc, 

c  x(5)=k2,  x(6)=k. 

c  Complete  Aerodynamics  (symmetric 

c  aeroforce  curves) 

c  Warren  Chen,  4-29-93 
cccccccccc 

integer  n, np, indx, i , t , j 
parameter  (nn=6) 

real*8  x(l:nn) ,fvec(l:nn) , f jac (1 :nn, 1 :nn) , delta (nn) 
$  , fs, fq,avib,avibd, alphaO,alphaOd,philt,phil 
real  y(nn,nn) , d, thetaO , thetaOd, q2s 
write (*,*) 'Enter  Root  Angle  (in  degrees) ' 
read ( * , * ) thetaOd 
thetaO=thetaOd*3 . 14159/180 . 0 
write(*,*)  Enter  q2s:' 
read (*, *) q2s 

write (*,*) 'Enter  Initial  value  of  qlO:' 
read(*,  *)x(l) 

write (*,*) 'Enter  Initial  value  of  q20:' 
read(*,  *)x(2)  ‘ 

write (*,*) 'Enter  Initial  value  of  qls:' 
read(*,  *)x(3) 

write (*,*) 'Enter  Initial  value  of  qlc:' 
read(*,*)x(4) 

write (*,*) 'Enter  Initial  value  of  k2 : ' 
read (*, *) x(5) 

write (*,*) 'Enter  Initial  value  of  k:' 

read (* , *) x (6) 

n=6 

np=6 

do  i=l,n 
delta  (i) =5 . 0 
enddo 
t=0 

open ( 5 , f ile= ' harm7dl .in', status= ' old ' ) 

open ( 9 , f ile= ' harm7dl . out ' , status= ' unknown ' ) 

write  (9, *)  '  ' 

write  ( 9 , * )  ' Harm7dl . out ' 

write (9, *) 'Complete  aerodynamics' 


write  (9, *)  '  ' 

write (9, *) ' INPUT' 

write ( 9 , 3 000 )' Root  angle:  ',theta0d, '  deg' 
write (9, *) '  ' 

write (9,5500) 'q2s ' , 'qlO ' , 'q20 ' , ' qls ' , ' qlc ' , ' k2 ' 
$  , 'k' , 'Fs (m/s) ' 

$  , 'w(Hz) ' , 'aO' , 'aV 

write (9,*) ' - 

10  continue 

if  (dabs (delta (6) ). le . 1 . Od-4  .and. 

$  dabs (delta (5) ). le . 1 . Od-4  .and. 

$  dabs (delta (4) ). le . 1 . Od-4  .and. 

$  dabs (delta (3) ). le. 1 . Od-4  .and. 

$  dabs (delta (2) ). le . 1 . Od-4  .and. 

$  dabs (delta (1) ) .le.l. Od-4) then 
goto  500 

endif 

t=t+l 

call  funcv (thetaO , q2s , alphaO , avib, phil , philt 
$  ,n,x,fvec) 

call  fdjac (q2s, thetaO , n, x, fvec, np, f jac) 
do  i=l,n 
do  j=l,n 
y(i, j)=o. 
enddo 
y (i, i) =1 . 
enddo 

call  ludcmp (f jac, n, np, indx, d) 
do  j=l,n 

call  lubksb (f jac, n,np, indx,y (1, j ) ) 
enddo 
do  j=l,n 
delta (j ) =0 . 
do  i=l,n 

delta (j ) =delta ( j ) +y ( j , i) *fvec (i) 
enddo 
enddo 

do  i=l,n 

X (i) =x (i) -delta (i) 
enddo 

if  (t .ge.25) then 


goto  500 
endif 
goto  10 
500  continue 

if  (t  .ge.  25) then 
write (*,*) 'Not  converged' 
goto  600 
endif 

fs=10.882/x(5) 

fq= (f s*x (6) /O . 07) /6 . 28 

avibd=avib*180 . 0/3 . 14159 

alpha0d=alpha0*180 . 0/3 . 14159 

write  (9, 5000)  q2s,x(l)  ,x(2),x(3),x(4),x(5) 

$  ,x(6) , fs, fq,alpha0d, avibd 
do  i=l,n 

X (i) =x (i) +0 . 0001 
enddo 
do  i=l,n 
delta(i)=5.0 
enddo 
t=0 

read (5, *) q2s 

if  (q2s  .ge.  999.0)  then 
goto  600 
endif 
goto  10 
600  continue 

write (9, *) '  ' 
write (9, *)' End  of  data' 
stop 
c 

c  formats 

1000  formate  ' ,  i2,  lx,  i2,  a,  fl7 . 9) 

2000  formate  ',fl7.9) 

3000  formate  ',a,f8.3,a) 

4000  formate  ',a,f6.4) 

4200  formate  '  ,a8,a8,a8,a8,a8,a8) 

4400  formate  '  ,  f  9 . 4  ,  f  9 . 4  ,  f  9 . 4  ,  f  9 . 4  ,  f  8 . 4  ,  f  7 . 4) 

5000  formate  '  ,  f4 . 3  ,  f  7 . 4  ,  f7 .4  ,  f  7 . 4  ,  f  7 . 4  ,  f  7 . 4  ,  f  7 .4  ,  f  8 . 3 
$  ,f8.3,f7.3,f7.3) 
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5100  formate  '  ,  i3  ,  f  8 . 4  ,  f  8 . 4 ,  f  8 . 4  ,  f  8 . 4) 

5500  formate  '  ,  a3  ,  a6 ,  a6 ,  a7 ,  a7,  a7  ,  a7 ,  a8  ,  a8  ,  a7 ,  a7 ) 

5600  format ( '  ' , a, a7 , a7 , a7 , a7) 

end 

c  Subroutine  to  find  the  Jacobian  matrix 

subroutine  fdjac (q2s, theta0,n,x, fvec,np,df ) 
integer  n,np,nmax 

real*8  df (np,np) , fvec (n) ,x(n) ,eps,f (40) ,alpha0,avib 
$  ,phil,philt 

parameter  (nmax=40 , eps=l . e-4 ) 
integer  i,j 
real  h, temp, thetaO 
do  j=l,n 
temp=x ( j ) 
h=eps*abs (temp) 
if (h.eq. 0 . ) h=eps 
X ( j ) =temp+h 
h=x ( j ) -temp 

call  funcv ( thetaO , q2s , alphaO , avib, phil , phi It , n, x, f ) 
X ( j ) =temp 
do  i=l,n 

df ( i , j ) = ( f ( i ) - fvec ( i ) ) /h 
enddo 
enddo 
c 

return 

end 

c  subroutine  to  find  flutter  velocity,  frequency 

subroutine  funcv (thetaO , q2s , alphaO , avib, phil 
$  ,philt,n,x,P) 
integer  n 

real  M,  la, Pi, rho, c, 1 , b, wh, wa, q2c, q2s, II, 12 , 13 , 14 , 15 
$  , thetaO , alphal , alphalt , sll , sl2 , sl3 , aol 
$  , laml, lam2, sml, sm2, sm3 
$  , aom,bl,b2,phia,phih 
real*8  x(n) ,P(n) , om, ra, CL20 , CM2 0 , alphaO 
$  , alphas , alphac , avib , mu , f , G , BoL , BIL 
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$  , B2L, B3L, B4L, BoM, BIM, B2M, B3M, B4M, pent , pent 2 
$  , phil , philt , dcLO , dcLsl , rlL, r2L, r3L 
$  , klL, k2L, rslL, rclL, ALsl 
$  ,BLcl,-dcMO,dcMsl,rlM,  r2M,  r3M,  klM 
$  ,k2M,rslM,rclM,AMsl,BMcl 


c  Variable  initialization 


11=1.0 

12=0.6779 

13=0.5 

14=0.783 

15=0.6366 

q2c=0 . 0 

M=0.158 

Ia=192.0  e  -6 

Pi=3 . 14159 

rho=1.23 

c=0.14 

1=0.559 

b=c/2 . 0 

wh=27.02 

wa=154.6 

phia=0 . 844605 

phih=l. 31538 

alphal=0.2094 

alphalt=-0 . 2094 

alpha0=theta0+0 . 5*phia*x (2) 

alphas=phia* (0 . 5*q2s+0 . 25*x (6) *q2c) 

$  +phih* (x (6) *x (4) ) 
alphac=phia* (0 . 5*q2c-0 . 25*x (6) *q2s) - 
$  phih*  (x(6)  *x(3)  ) 
avib=dsqrt (alphas**2+alphac**2) 
om=wh/wa 

tnu=  (M/1)  /  (Pi*rho*b**2) 
ra= (dsqrt (la/M) ) /b 

c 

c  Linear  Lift 
sll=3 .142 
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Sl2=1.571 
sl3=0.0 
aol=5.73 
laml=0  rl5 
lam2=0.55 

F= (x(6) **2*lam2+laml**2) / (laml**2+x (6 )  **2) 

G= (x (6) *laml* (lam2-l) ) / (laml**2+x (6) **2) 

BoL=x (2 ) *aol*0 . 5 

BlL=x(6) * (0.25*x(6) *sll-0 . 5*x (6) *sl2 
$  -0.25*G*aol)+0.5*F*aol 
B2L=x(6) * (0 . 5*sll+0 . 5*sl3+0 . 25*F*aol) +0.5*G*aol 
B3L=x(6) * (sll*x(6) +G*aol) 

B4L=F*aol*x(6) 

c  Linear  Moment 

sml=-0 . 786 
sm2=-0 . 589 
sm3=-0 . 786 
aom=0 . 0 

BoM=x ( 2 ) *aom*  0 . 5 

BlM=x (6) * (0 . 25*x (6) *sMl-0 . 5*x (6) *sM2-0 . 25*G*aoM) 
$  +0.5*F*aoM 

B2M=x(6) * (0 . 5*sMl+0 . 5*sM3+0 .25*F*aoM) +0 . 5*G*aoM 
B3M=x(6) * (sMl*x(6) +G*aoM) 

B4M=F*aoM*x(6) 

c 

c  Nonlinear  Lift 
bl=8.6 

pent= (alphal-alphaO) /avib 
if  (pent  .gt.  1.0)  then 
phil=1.571 

else  if  (pent  .It.  -1.0)  then 

phil=-1.571 

else 

phil=dasin (pent) 
endif 
c 

pent2= (alphalt-alphaO) /avib 
if  (pent 2  .gt.  1.0)  then 
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philt=1.571 

else  if  (pent2  .It.  -l.O)  then 

philt=-1.571 

else 

philt=dasin (pent2 ) 
endif 
c 

dcLO= (bl*avib/3 . 14159) * ( -pent* (1 . 571- 
$  phil) +dcos (phil) ) - (bl*avib/3 . 14159) * 

$  (pent2* (philt+1 . 571) +dcos (philt) ) 
dcLsl= (bl*avib/3 . 14159) * ( (1.571-phil) - 
$  0 . 5*dsin (2 . 0*phil) ) 

$  - (bl*avib/3 . 14159) * {- (philt+l . 571) - 
$  0 . 5*dsin (2 . 0*philt) ) 

CL20=-dcL0 

rlL=0 . 6+0 . l*dcL0**2 

r2L= (0.2465  +  0 . 005*dcL0**2 )  **2 

r3L= (-3 . 7399+7 . 0*dcL0**2) *r2L 

]clL=r2L-x(6)  *x(6) 

]c2L=rlL*x(6) 
rslL=-r2L*dcLsl 
rclL=-r3L*x (6) *dcLsl 

C 

ALsl= {klL*rslL+k2L*rclL) / (klL**2+k2L**2) 
BLcl= {klL*rclL-k2L*rslL) / (klL**2+k2L**2) 
c 
c 

c  Nonlinear  Moment 
b2=0.48 

dcM0= (b2*avib/3 .14159) * (-pent* (1.571- 
$  phil) +dcos (phil) ) 

$  - (b2*avib/3 .14159) 

$  * (pent2* (philt+1 . 571) +dcos (philt) ) 
dcMsl=(b2*avib/3. 14159) * ( (1.571-phil) 

$  -0 . 5*dsin (2 . 0*phil) ) 

$  - (b2*avib/3 .14159) * (- (philt+1 . 571) - 
$  0 . 5*dsin (2 . 0*philt) ) 

CM20=-dcM0 

rlM=0 . 25+0 . l*dcM0**2 
r2M= (2 . 0+0 . l*dcM0**2) **2 
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r3M= (5.0-0 . 6*dcM0**2) *r2M 
klM=r2M-x{6) *x(6) 
k2M=rlM*x(6) 
rslM=-r2M*dcMsl 
rclM=-r3M*x (6) *dcMsl 
c 

AMsl= (klM*rslM+k2M*rclM) / (klM**2+k2M**2 ) 

BMcl= (klM*rclM-k2M*rslM) / (klM**2+k2M**2 ) 

c  Equations 

p (1) =mu*Pi*Il*om**2*x (5) **2*x(l) -I4*aol*theta0 
$  -I2*BoL-I4*CL20 
c 

p{2)=mu* (Pi/4.0) *ra**2*I3*x(5) **2*x(2) -0 . 25*I5*aol* 
$  thetaO-0 .25*I3*BoL-0 . 25*I5*CL20-I5*CM20 
c 

p (3) =mu*Pi*Il* (-x(6) **2*x(3) +om**2*x(5) **2*x(3) ) 

$  -12* (q2s*BlL-q2c*B2L) -II* (x(3) *B3L+x(4) *B4L) 

$  -l4*ALsl*alphas/avib+I4*BLcl*alphac/avib 
c 

p (4) =mu*Pi*Il* (-X (6) **2*x (4) +om**2*x(5) **2*x(4) ) 

$  -12* (q2s*B2L+q2c*BlL) -II* (-x(3) *B4L+x(4) *B3L) 

$  -I4*ALsl*alphac/avib-I4*BLcl*alphas/avib 

c 

p (5) =mu* (Pi/4.0) *ra**2*I3* (-x(6) **2*q2s 
$  +x(5) **2*q2s) 

$  -0.25*13* (q2s*BlL-q2c*B2L) -0.25*12* (x(3) *B3L 
$  +x(4)*B4L) 

$  -13* (q2s*BlM-q2c*B2M) -12* (x(3) *B3M+x (4) *B4M) 

$  -I5*AMsl*alphas/avib+I5*BMcl*alphac/avib 
$  -0 .25*I5*ALsl*alphas/avib+0 . 25*I5*BLcl*alphac/avib 
c 

p (6) =mu* (Pi/4 .0) *ra**2*I3* (-X (6) **2*q2c 
$  +x(5) **2*q2c) 

$  -0.25*13* (q2s*B2L+q2c*BlL) -0 .25*12* (-x(3) *B4L 
$  +x(4)*B3L) 

$  -13* (q2s*B2M+q2c*BlM) -12* (-x(3) *B4M+x(4) *B3M) 

$  -I5*AMsl*alphac/avib-I5*BMcl*alphas/avib 
$  -0 .25*I5*ALsl*alphac/avib-0 .25*I5*BLcl*alphas/avib 


return 
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end 

c 

c 

c 

c  L-U  decomposition 

subroutine  ludcmp (a, n, np, indx, d) 
parameter  (nmax=100 , tiny=l . Oe-20) 
double  precision  a (np, np) , w (nmax) 
dimension  indx(n) 

d=l. 

do  12  i=l,n 
aamax=0 . 
do  11  j=l,n 

if  (dabs  (a  (i, j ) )  .gt .aamax)  aamax=dabs (a (i , j ) ) 

11  continun 

if  (aamax. eq. 0 . )  pause  'singular  matrix.' 
w ( i ) =1 . /aamax 

12  continue 

do  19  j=l,n 
do  14  i=l,j-l 
sum=a (i , j ) 

do  13  k=l,i-l 

sum=sum-a (i,k)*a(k,j) 

13  continue 
a (i, j ) =sum 

14  continue 
aamax=0 . 

do  16  i=j,n 
sum=a (i, j ) 
do  15  k=l,j-l 

sum=sum-a (i, k) *a  (k, j ) 

15  continue 
a (i, j ) =sum 

dum=vv(i) *abs (sum) 
if  (dum.ge. aamax)  then 
imax=I 
aamax=dum 
endif 
continue 

if  (j .ne.imax) then 
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do  17  k=l,n 
dum=a (imax, k) 
a (imax, k) =a ( j , k) 
a  ( j  ,~k)  =dum 
continue 
d=-d 

w  (imax)  =w  ( j  ). 
endif 

indx( j ) =imax 

if(a(j,j) .eq.0.)a(j,j) =tiny 
if ( j .ne .n) then 
dum=l . /a ( j , j ) 
do  18  i=j+l,n 
a(i,j)=a(i,j) *dum 
continue 
endif 
continue 
return 
end 


Back  substitution 

subroutine  lubksb (a, n, np, indx, b) 
dimension  indx(n),b(n) 
double  precision  a(np,np) 
ii=0 

do  12  i=l,n 
ll=indx (i) 
sum=b (11) 
b(ll)=b(i) 
if  (ii .ne . 0) then 
do  11  j=ii,i-l 

sum=sum-a (i, j) *b(j) 
continue 

else  if  (sum.ne.O.)  then 
ii=I 
endif 
b(i)=sum 
continue 
do  14  i=n, 1, -1 
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sum=b ( i ) 
do  13  j=i+l,n 
sum=sum-a (i, j ) *b ( j ) 
continue 
b (i) =sum/a (i, i) 
continue 
return 
end 
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cccccccccc 

c  staticsl  .-FOR  Statics  divergence  program, 
c  x(l)=qlO,  x(2)=q20 

c  Complete  Aerodynamics 

c  A  built-in  loop  to  read  k2 ' s  from  staticsl. in, 

c  perform  Newton  Raphson  Method  and  takes  the 

c  resultant  qlO  and  q20  as  the  initial  guess  for 

c  the  next  velocity  to  insure  convergence 

c  to  the  correct  roots, 

c  Warren  Chen,  4-29-93 
cccccccccc 

integer  n, np, indx, i , t , j 
parameter  (nn=2) 

real *8  x(l:nn) ,fvec(l:nn) , f jac (1 :nn, 1 :nn) , delta (nn) 
$  ,alphal,aoad,aoa,v 
real  y(nn,nn) , d, thetaO , thetaOd, kat 
write (*,*)' Enter  Root  Angle  of  attack' 
read (* , * ) thetaO 

write (*,*) 'Enter  the  initial  guess  for  qlO:' 
read(*,  *)x(l) 

write (*,*) 'Enter  the  initial  guess  for  q20:' 
read(*,  *)x(2) 
alphal=0 . 2094 

theta0d=theta0*180 . 0/3 . 14159 

n=2 

np=2 

do  i=l,n 
delta (i) =5 . 0 
enddo 
t=0 

open (5, f ile= ' staticsl . in' , status= ' old ' ) 
open(9, file= ' staticsl. out ' , status=' unknown' 

$  ,  access= ' sequential ' ) 
read (5, *)kat 
write (9,*)'  ' 

write ( 9 , * ) ' staticsl . out ' 
write (9 ,*)' Complete  aerodynamics' 
write (9 , ' 

write (9, 5600) 'Theta-R' , 'k2' , 'qlO' , 'q20' 

$  , 'AOA{deg) ' , 'V(m/s) ' 


I 


write (9,*)' - 

10  continue 

if  (dabs (delta (2) ). le. 1 . Od-3  .and. 

$  dabs (delta (1) ). le . 1 . Od-3) then 
goto  500 
endif 
t=t+l 

call  funcv (kat , thetaO , n, X, fvec) 
call  f djac (kat, thetaO , n, x, fvec, np, f jac) 
do  i=l, n 
do  j=l,n 
y(i, j)=0. 
enddo 
y(i,i)=l. 
enddo 

call  ludcmp (f jac, n, np, indx, d) 
do  j=l,n 

call  lubksb (f jac, n, np, indx,y (1, j ) ) 
enddo 
do  j=l,n 
delta ( j ) =0 . 
do  i=l,n 

delta (j ) =delta ( j ) +y ( j , i) *fvec (i) 
enddo 
enddo 

do  i=l,n 

X (i) =x(i) -delta (i) 
enddo 

if  (t.ge.25) then 
goto  500 
endif 
goto  10 
500  continue 

aoa=theta0+0 . 5*0 . 8446*x (2) 
aoad=aoa*180/3 . 14159 
v=154.6*0.07/kat 

write (9, 5100) thetaO, kat, x(l) ,x(2) ,aoad,v 

x(l)=x(l)+0.0001 

x{2) =x(2) +0.0001 

delta (1) =5 . 0 

delta (2) =5.0 


t  =  0 

read (5 , *) kat 

if  (kat  .ge,  999.0)  then 
goto  600 
endif 
goto  10 

600  continue 

write (9 , *)  '  ' 

write (9, *)' End  of  Data' 

stop 

c  formats 

1000  format ( '  ' , i2 , lx, i2 , a, f 17 . 9) 

2000  format ( '  ',fl7.9) 

3000  format ( '  ',a,f8.3,a) 

4000  format ( '  ',a,f6.3) 

4200  f ormat ( '  ' , a8 , a8 , a8 , a8 , a8 , a8) 

4400  formate  '  ,  f  9 . 4  ,  f  9 . 4 ,  f  9 . 4  ,  f  9 . 4  ,  f  8 . 4  ,  f  7 . 4 ) 
5000  formate  '  ,  i3  ,  f  9 . 4  ,  f  9 . 4  ,  f  9 . 4  ,  f  9 . 4  ,  f  8 . 4  ,  f  7 . 4 ) 
5100  formate  '  ,  f  7 . 3  ,  f  7 . 3  ,  f  7 . 4 ,  f  7 . 4  ,  f  9 . 3  ,  f  9 . 3 ) 
5500  f ormat ( '  ' , a, a8 , a8 , a8 , a8 , a8 , a8) 

5600  f ormat  ( '  ' , a7 , a7 , a7 , a7 , a9 , a9 ) 

end 


subroutine  fdjac (kat , theta0,n,x, fvec,np,df ) 
integer  n,np,nmax 

real *8  df (np, np) , fvec (n) , x (n) , eps, f (40) 
parameter  (nmax=40 , eps=l . e-4 ) 
integer  i,j 

real  h, temp, kat , thetaO 
do  j=l,n 
temp=x ( j ) 
h=eps*abs (temp) 
if (h.eq. 0 . ) h=eps 
x ( j ) =temp+h 
h=x ( j ) - temp 

call  funcv(kat,thetaO,n,x,f) 
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X ( j ) =temp 
do  i=l,n 
df(i,j)  =  (f(i: 
enddo 
enddo 
return 
end 


•fvec (i) ) /h 


subroutine  funcv (kat , thetaO , n, x, P) 
integer  n 

real  M,  la, Pi , rho, 1 , wh, wa, II, 13 , 14 , 15 , kat 
$  ,bl,b2,phia,phih, 

$  alphaO,alphal,alphav, thetaO,b,c, 12, mu 
real*8  x(n) ,P(n) ,om, ra,CL20,CM20,aol,pent,BoL 
aol=5.73 
alphal=0 .2094 
alphav=0 .002 
11=1.0 
12=0.6779 
13=0.5 
14=0.783 
15=0 . 6366 
M=0.158 
Ia=192.0  e  -6 
Pi=3 .14159 
rho=l .23 
c=0.14 
1=0.559 
b=c/2.0 
wh=27. 02 
wa=154.6 
om=wh/wa 

mu= (M/1) / {Pi*rho*b**2) 
ra= (dsqrt (la/M) ) /b 
BoL=x (2) *aol*0 . 5 

c  nonlinear  lift 

bl=8.6 
phia=0 . 8446 
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phih=l .3154 

alpha0=theta0+0 . 5*phia*x (2) 
pent= (alphal-alphaO) /alphav 
if  (pent  .gt.  1.0)  then 
phil=1.571 

else  if  (pent  .It.  -1.0)  then 

phil=-l . 571 

else 

phil=dasin (pent) 
endif 

CL20=- (bl*alphav/3 .14159) * (-pent* (1.571- 
$  phil) +dcos (phil) ) 

c  nonlinear  moment 

b2=0 .48 

CM20=- (b2*alphav/3 . 14159) * ( -pent* (1.571- 
$  phil) +dcos (phil) ) 

c  equations 

p (1) =mu*Pi*Il*om**2*kat**2*x (1) -I4*aol*theta0 
$  -I2*BoL-I4*CL20 
c 

p (2) =mu* (Pi/4 . 0) *ra**2*I3*kat**2*x(2) - 
$  0 .25*I5*aol*theta0 

$  -0 . 25*I3*BoL-0 . 25*I5*CL20-I5*CM20 

c 

return 

end 


c  L-U  decomposition 

subroutine  ludcmp (a,n,np, indx,d) 
parameter  (nmax=100 , tiny=l . Oe-20) 
double  precision  a (np, np) , w (nmax) 
dimension  indx(n) 
d=l. 

do  12  i=l,n 
aamax^O . 
do  11  j=l,n 


if  (dabs (a (i, j ) ) .gt . aamax)  aamax=dabs (a (i, j ) ) 
continun 

if  (aamax. eq. 0 . )  pause  'singular  matrix.' 
w(i)  =1.  /aamax 
continue 
do  19  j=l,n 
do  14  i=l,j-l 
sum=a (i, j ) 

do  13  k=l,i-l 

sum=sum-a {i,k) *a (k,  j  ) 
continue 
a (i , j ) =sum 
continue 
aamax=0 . 
do  16  i=j,n 
sum=a (i, j ) 
do  15  k=l,j-l 

sum=sum-a (i , k) *a (k,  j  ) 
continue 
a (i, j ) =sum 
dum=w  ( i )  *abs  ( sum) 
if  (dum.ge. aamax)  then 
imax=I 
aamax=dum 
endif 
continue 

if  ( j . ne . imax) then 
do  17  k=l,n 
dum=a (imax, k) 
a (imax,k) =a ( j ,k) 
s ( j , k) =dum 
continue 
d=-d 

w(imax)  =w(j) 
endif 

indx(j)  =imax 

if(a(j,j) .eq.0.)a(j,j) =tiny 
if (j .ne.n) then 
dum=l./a(j, j) 
do  18  i=j+l,n 
a(i,j)=a(i,j) *dum 
continue 
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endif 

continue 

return 

end 


c  Back  substitution 

subroutine  lubksb (a, n, np, indx,b) 
dimension  indx(n),b(n) 
double  precision  a(np,np) 
ii=0 

do  12  i=l,n 
ll=indx (i) 
sum=b (11) 
b(ll)=b(i) 
if  (ii .ne . 0) then 
do  11  j=ii,i-l 

sum=sum-a (i , j ) *b ( j ) 

11  continue 

else  if  (sum.ne.O.)  then 
ii=I 
endif 
b (i) =sum 

12  continue 

do  14  i=n, 1,-1 
sum=b(i) 
do  13  j=i+l,n 
sum=sum-a (i, j ) *b ( j ) 

13  continue 
b (i) =sum/a(i, i) 
continue 
return 
end 
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